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First,  we  describe  an  e-DFT  protocol  in  which  the  non-additive  kinetic  energy 
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implementation  of  the  exact  calculation  of  the  non-additive  kinetic  potential  (NAKP) 
and  apply  it  to  molecular  systems.  We  demonstrate  that  the  implementation  using 
the  exact  NAKP  is  in  excellent  agreement  with  reference  Kohn-Sham  calculations, 
whereas  the  approximate  functionals  lead  to  qualitative  failures  in  the  calculated 
energies  and  equilibrium  structures. 

Next,  we  introduce  density-embedding  techniques  to  enable  the  accurate  and  stable 
calculation  of  correlated  wavefunction  (CW)  in  complex  environments.  Embedding 
potentials  calculated  using  e-DFT  introduce  the  efcct  of  the  environment  on 
a  subsystem  for  CW  calculations  (WFT-in-DFT).  We  demonstrate  that  WFT-in- 
DFT  calculations  are  in  good  agreement  with  CW  calculations  performed  on  the  full 
complex. 
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between  subsystems  by  introduction  of  a  projection  operator.  Utilizing  the  projection 
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Methods  that  exploit  the  intrinsic  locality  of  molecular  interactions  show  significant 
promise  in  making  tractable  the  electronic  structure  calculation  of  large-scale  sys¬ 
tems.  In  particular,  embedded  density  functional  theory  (e-DFT)  offers  a  formally 
exact  approach  to  electronic  structure  calculations  in  which  the  interactions  between 
subsystems  are  evaluated  in  terms  of  their  electronic  density.  In  the  following  disser¬ 
tation,  methodological  advances  of  embedded  density  functional  theory  are  described, 
numerically  tested,  and  applied  to  real  chemical  systems. 

First,  we  describe  an  e-DFT  protocol  in  which  the  non-additive  kinetic  energy 
component  of  the  embedding  potential  is  treated  exactly.  Then,  we  present  a  general 
implementation  of  the  exact  calculation  of  the  non-additive  kinetic  potential  (NAKP) 
and  apply  it  to  molecular  systems.  We  demonstrate  that  the  implementation  using 
the  exact  NAKP  is  in  excellent  agreement  with  reference  Kohn-Sham  calculations, 
whereas  the  approximate  functionals  lead  to  qualitative  failures  in  the  calculated 
energies  and  equilibrium  structures. 

Next,  we  introduce  density-embedding  techniques  to  enable  the  accurate  and  sta¬ 
ble  calculation  of  correlated  wavefunction  (CW)  in  complex  environments.  Embed¬ 
ding  potentials  calculated  using  e-DFT  introduce  the  effect  of  the  environment  on 
a  subsystem  for  CW  calculations  (WFT-in-DFT).  We  demonstrate  that  WFT-in- 
DFT  calculations  are  in  good  agreement  with  CW  calculations  performed  on  the  full 
complex. 

We  significantly  improve  the  numerics  of  the  algorithm  by  enforcing  orthogonality 
between  subsystems  by  introduction  of  a  projection  operator.  Utilizing  the  projection- 


based  embedding  scheme,  we  rigorously  analyze  the  sources  of  error  in  quantum 
embedding  calculations  in  which  an  active  subsystem  is  treated  using  CWs,  and  the 
remainder  using  density  functional  theory.  We  show  that  the  embedding  potential 
felt  by  the  electrons  in  the  active  subsystem  makes  only  a  small  contribution  to  the 
error  of  the  method,  whereas  the  error  in  the  nonadditive  exchange-correlation  energy 
dominates.  We  develop  an  algorithm  which  corrects  this  term  and  demonstrate  the 
accuracy  of  this  corrected  embedding  scheme. 
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Summary 

Quantum-chemistry  calculations  invariably  feature  a  compromise  between  computa¬ 
tional  efficiency  and  accuracy.  Kohn-Sham  Density  Functional  Theory  (KS-DFT)  is 
one  of  the  most  commonly  used  methods  because  it  is  tractable  for  systems  contain¬ 
ing  up  to  a  thousand  atoms.  However,  DFT  methods,  which  generally  employ  ap¬ 
proximate  descriptions  of  electron  correlation,  often  fail  to  even  qualitatively  predict 
reaction  barriers,  which  are  critical  in  determining  reaction  rates  and  mechanisms. 
More  rigorous  correlated  wavefunction  (CW)  methods  can  provide  an  accurate  and 
systematically  improvable  description  of  reaction  barriers,  but  poor  scaling  of  these 
methods  can  lead  to  impractical  computation  costs  for  systems  with  even  tens  of 
atoms. 

Embedded  density  functional  theory  (e-DFT),  is  a  formally  exact  approach  to 
electronic  structure,  which  exploits  the  locality  of  molecular  interactions.  In  the  e- 
DFT  approach,  a  system  is  divided  into  smaller  subsystems.  Then,  the  interactions 
between  subsystems  are  evaluated  in  terms  of  their  electronic  density.  This  approach 
naturally  leads  to  two  different  partitioning  schemes.  One,  a  large  system  can  be  sep¬ 
arated  into  many  smaller  subsystems,  allowing  for  significant  computational  savings 
as  this  leads  to  a  naturally  parallclizable  algorithm.  Two,  a  large  system  is  sepa¬ 
rated  into  a  ‘high-level’  (accurate)  and  a  ‘low-level’  (approximate)  regions;  therefore, 
if  higher  accuracy  is  required  for  a  region  of  system,  a  higher  level  of  theory  can 
be  seamlessly  embedded  into  a  DFT  environment  and  still  remain  computationally 
tractable  as  only  a  handful  of  atoms  are  treated  at  the  CW  level.  The  objectives  of  the 
e-DFT  approach  are  thus  similar  to  those  of  more  approximate  partitioning  and  frag- 
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mentation  schemes;  however,  e-DFT  avoids  the  uncontrolled  approximations  (such  as 
link  atoms)  and  errors  associated  with  subsystem  interfaces  that  fundamentally  limit 
other  widely  used  methods. 

This  dissertation  is  focused  on  the  development  of  new  methods  that  make  e-DFT 
an  accurate,  practical,  and  scalable  method  for  the  description  of  complex  systems. 
The  following  key  contributions  are  discussed:  (i)  the  development  of  numerically 
exact  methods  for  obtaining  subsystem  embedding  potentials  in  e-DFT,  which  reduce 
embedding  errors  by  orders  of  magnitude  in  comparison  with  previously  available 
approximate  methods,  (ii)  the  development  of  parallelization  algorithms  that  enable 
the  description  of  large  systems  with  sub-linear  scaling  of  the  required  computational 
time,  (iii)  the  combination  of  e-DFT  potentials  with  correlated  wavefunction  theory 
(WFT)  methods  to  enable  seamless  WFT-in-DFT  embedding  for  general  systems, 
and  (iv)  the  development  of  systematically  improvable  methods  for  WFT-in-DFT  to 
allow  for  a  hierarchy  of  embedding  methods  with  increasing  accuracy. 

In  Chapter  1,  we  describe  an  embedded  density  functional  theory  (DFT)  proto¬ 
col  in  which  the  non-additive  kinetic  energy  component  of  the  embedding  potential 
is  treated  exactly.  At  each  iteration  of  the  Kohn-Sham  equations  for  constrained 
electron  density,  the  Zhao-Morrison-Parr  constrained  search  method  for  constructing 
Kohn-Sham  orbitals  is  combined  with  the  King-Handy  expression  for  the  exact  kinetic 
potential.  We  use  this  exact  embedding  protocol  to  calculate  ionization  energies  for 
a  series  of  3-  and  4-electron  systems,  and  the  results  are  compared  to  embedded  DFT 
calculations  that  utilize  the  Thomas- Fermi  (TF)  and  the  Thomas-Fcrmi-von  Wei- 
sacker  approximations  to  the  kinetic  energy  functional.  These  calculations  illustrate 
a  breakdown  clue  to  the  TF  approximation  for  the  non-additive  kinetic  potential,  with 
errors  of  30-80%  in  the  calculated  ionization  energies;  by  contrast,  the  exact  protocol 
is  found  to  be  accurate  and  stable.  This  work  has  been  published  as  J.  D.  Goodpaster, 
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N.  Ananth,  F.  R.  Manby,  and  T.  F.  Miller  III,  “Exact  non-additive  kinetic  potentials 
for  embedded  density  functional  theory,”  J.  Chem.  Phys.,  133,  084103  (2010). 

In  Chapter  2,  we  present  a  general  implementation  of  the  exact  calculation  of 
the  non-additive  kinetic  potential  (NAKP).  Potential  energy  curves  are  computed 
for  the  dissociation  of  Li+-Be,  CH3-CF3,  and  hydrogen-bonded  water  clusters,  and 
e-DFT  results  obtained  using  this  exact  treatment  of  the  NAKP  are  compared  with 
those  obtained  using  approximate  kinetic  energy  functionals.  In  all  cases,  the  exact 
NAKP  is  in  excellent  agreement  with  reference  Kohn-Sham  calculations,  whereas  the 
approximate  functionals  lead  to  qualitative  failures  in  the  calculated  energies  and 
equilibrium  structures.  We  also  demonstrate  an  accurate  pairwise  approximation  to 
the  NAKP  that  allows  for  efficient  parallelization  of  the  method  in  large  systems; 
benchmark  calculations  on  molecular  crystals  reveal  ideal,  size-independent  scaling 
of  wall-clock  time  with  increasing  system  size.  This  work  has  been  published  as  J. 
D.  Goodpaster,  T.  A.  Barnes,  and  T.  F.  Miller  III,  “Embedded  density  functional 
theory  for  covalently  bonded  and  strongly  interacting  subsystems,”  J.  Chem.  Phys., 
134,  164108  (2011). 

In  Chapter  3,  we  introduce  density  embedding  techniques  to  enable  the  accu¬ 
rate  and  stable  calculation  of  CWs  in  complex  environments.  Embedding  potentials 
calculated  using  e-DFT  introduce  the  effect  of  the  environment  on  a  subsystem  for 
wavefunction  calculations  (WFT-in-DFT).  These  methods  are  demonstrated  for  cal¬ 
culating  the  potential  energy  curve  of  the  dispersion-bound  ethene-propene  dimer 
and  for  the  hexaaquairon(II)  ion  (HA).  The  potential  energy  curve  for  the  ethene- 
propene  dimer  reveals  that  benchmark  CCSD(T)  calculations  performed  over  the  full 
complex  can  be  reproduced  within  0.05  kcal/mol  using  this  method,  illustrating  that 
small  energy  differences  can  be  accurately  calculated  while  embedding  across  a  cova¬ 
lent  bond.  e-DFT  calculations  on  HA  that  employ  these  new  techniques  demonstrate 
that  Kohn-Sham  (KS-DFT)  calculations  can  be  reproduced  for  both  the  low-spin  and 
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the  high-spin  states.  The  ability  of  different  exchange-correlation  (XC)  functionals  to 
describe  the  energy  differences  between  the  low-spin  and  the  high-spin  states  (AFAh) 
and  to  describe  the  ligation  energy  in  HA  is  studied.  KS-DFT  calculations  of  AFAh 
demonstrate  a  strong  dependency  on  XC  functionals  where  WFT-in-DFT  calcula¬ 
tions  reveal  a  significantly  diminished  dependency  and  that  AFAh  is  in  good  agree¬ 
ment  with  ab  initio  calculations  performed  on  the  full  complex.  The  ligation  energies 
have  a  small  dependency  on  the  XC  functional  and  are  near  identical  for  KS-DFT 
and  WFT-in-DFT  demonstrating  that  the  interactions  energies  between  subsystems 
remain  at  DFT  level  accuracy.  This  work  has  been  published  as  J.  D.  Goodpaster,  T. 
A.  Barnes,  F.  R.  Manby,  and  T.  F.  Miller  III,  “Density  functional  theory  embedding 
for  correlated  wavefunctions:  Improved  methods  for  open-shell  systems  and  transition 
metal  complexes,”  J.  Chem.  Phys.,  137,  224113  (2012). 

In  Chapter  4,  we  present  continuing  work  on  projection  based  embedding.  The 
original  formulation  has  been  published  as  F.  R.  Manby,  M.  Stella,  J.  D.  Goodpaster, 
and  T.  F.  Miller  III,  “A  simple,  exact  density-functional-theory  embedding  scheme,” 
J.  Chem.  Theory  Comput.,  8,  2564  (2012)  and  T.  A.  Barnes,  J.  D.  Goodpaster, 
F.  R.  Manby,  and  T.  F.  Miller  III,  “Accurate  basis  set  truncation  for  wavefunction 
embedding,”  J.  Chem.  Phys.,  139,  024103  (2013).  Building  upon  that  work,  we 
analyze  the  sources  of  error  in  quantum  embedding  calculations  in  which  an  active 
subsystem  is  treated  using  CW  methods,  and  the  remainder  using  DFT.  We  show 
that  the  embedding  potential  felt  by  the  electrons  in  the  active  subsystem  makes  only 
a  small  contribution  to  the  error  of  the  method  whereas  the  error  in  the  nonadditive 
exchange-correlation  energy  dominates.  We  test  an  MP2  correction  for  this  term 
and  demonstrate  that  the  corrected  embedding  scheme  accurately  reproduces  the 
full  wavefunction  calculations  for  a  series  of  chemical  reactions.  This  work  has  been 
published  as  J.  D.  Goodpaster,  T.  A.  Barnes,  F.  R.  Manby,  and  T.  F.  Miller  III, 
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“Accurate  and  systematically  improvable  density  functional  theory  embedding  for 
correlated  wavefunctions,”  J.  Chem.  Phys.,  140,  18A507  (2014). 
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Chapter  1 

Exact  non-additive  kinetic  potentials  for  embedded 
density  functional  theory 

1.1  Introduction 

Orbital-free  embedded  density  functional  theory  (e-DFT)  is  an  appealing  method 
for  calculating  the  electronic  structure  of  complex  molecular  systems.  It  provides  a 
formally  exact  framework  for  dividing  the  total  electronic  density  of  a  system  into 
subsystem  densities  that  can  be  separately  calculated. 1  1  This  feature  of  e-DFT  allows 
for  the  development  of  multiscale  strategies  in  which  the  electronic  density  for  regions 
of  central  interest  is  calculated  with  high  accuracy,  while  the  electronic  density  for 
surrounding  regions  is  obtained  using  more  approximate  techniques.5- 

However,  in  addition  to  the  usual  approximations  for  the  basis  set  and  the  exchange- 
correlation  functional  that  appear  in  Kohn-Sham  (KS)  DFT,8  e-DFT  requires  the 
evaluation  of  a  non-additive  contribution  to  the  kinetic  energy  from  the  subsystem 
densities.4  This  term,  which  is  generally  largest  for  cases  in  which  the  subsystem 
densities  are  strongly  overlapping,9  is  a  significant  source  of  error  in  many  e-DFT 
calculations,  and  it  currently  limits  the  method  to  applications  in  which  the  subsys¬ 
tem  densities  involve  non-bonded  or  weakly  interacting  molecular  groups. 4,9  Although 
encouraging  progress  towards  the  accurate  calculation  of  the  non-additive  kinetic  en¬ 
ergy  contribution  have  been  reported, 3,9-13  more  work  in  this  direction  is  needed. 

In  this  paper,  we  present  an  exact  protocol  for  calculating  the  non-additive  kinetic 
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energy  contribution  in  e-DFT  calculations,  and  we  report  calculations  in  which  the 
protocol  is  applied  to  3-  and  4-electron  systems  that  exhibit  strongly  overlapping 
subsystem  densities.  Although,  in  its  numerically  demonstrated  form,  this  protocol 
does  not  offer  any  computational  speed-up  over  the  KS  calculation  for  the  full  system, 
it  suggests  new  methods  to  systematically,  efficiently,  and  accurately  perform  e-DFT 
calculations  for  large  systems,  which  we  discuss. 

1.2  Orbital- Free  Embedded  DFT 

Suppose  that  the  entire  electronic  density  pab  for  a  closed-shell  system  is  divided 
into  two  subsystems,  pa  and  pB,  such  that  pab  =  Pa  +  Pb-  The  one-electron  orbitals 
that  give  rise  to  these  subsystem  electronic  densities  obey  the  coupled  Kohn-Sham 
equations  for  constrained  electron  density  (KSCED), 


-iv2  +  KfCE>A,PB;r] 


AiA, 


[r)  =  ei  <p. 


i  =  (!■!) 


--V^  +  KrnPB,pA;r] 


»?(r)  =  4V?  W, 


Nb 

'  1 .  •  (1-2) 


where  N\  and  Nb  are  the  number  of  electrons  in  the  respective  subsystems, 


Na/2 

Pa(t)  =  2  |0?A(r)|2,  and 

i=  1 

Nb/ 2 

pB(r)  =  2  ^  |0f(r)|2. 


(1.3) 

(1.4) 


i= 1 

In  these  coupled  equations,  V^sced[pa,  Pb;  r]  is  the  KS  effective  potential  for  subsys¬ 
tem  A  embedded  in  subsystem  B, 


KffSCED[PA,  Pb;  r]  =  vne(r)  +  Vj[pAB ;  r]  +  vxc[pAB ;  r]  +  vnad[pA,  pB;  r], 


(1.5) 


and  V^gSCED[pB,  Pa;  r]  is  the  similarly  defined  KS  effective  potential  for  subsystem  B 
embedded  in  subsystem  A.  The  contributions  to  the  KS  effective  potential  include 


Vnuc 


1  v  nuc  ry 

”"’(r)  =  £|r-Ri|  * 

?.  1  ! 

(1.6) 

vj[pAB;r]  =  /  ,  ,  .dr,  and 

(1.7) 

J  |r  -  r| 

r  n  5EXC  [p] 

vxc[Pab]t\  = 

Op  P=PAB_ 

(r), 

(1.8) 

which  are  the  usual  nuclear-electron  Coulomb  potential,  Hartree  potential,  and  exchange- 
correlation  potential,  respectively,  and  NnvLC  is  the  number  of  nuclei  in  the  system. 
The  final  term  in  V^esced[pb,  Pa;  r]  is  the  non-additive  kinetic  potential  (NAKP) 


Aiad  [pa  5  Pb  1 1". 


'hrsnad[pA,pB]' 

M  _  6T,\p] 

<v\ 

&T,[p] 

dpA 

'  '  5p 

A  ) 

P=PAB 

5p 

(1.9) 


which  is  obtained  from  the  functional  derivative  of  the  non-additive  component  of  the 
non-interacting  kinetic  energy 


rsnad[pA, Pb]  =  Ts[pAB ]  -  Ts[pA]  -  Ts[pB}.  (1.10) 


The  total  energy  functional  for  the  embedded  system  is 


E[pAB\  =  Ts[pA]  +  Ts  [pB]  +  Tsliad[pA,pB]  +  Ke  [Pab]  +  J[p  ab]  +  Exc[pAB],  (1-11) 
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where  the  last  three  terms  on  the  right  hand  side  (RHS)  are  the  nuclear-electron 
Coulomb  energy,  Hartree  energy,  and  exchange-correlation  energy  for  the  total  den¬ 
sity. 

Two  aspects  of  the  orbital-free  embedding  DFT  formulation  are  worth  emphasiz¬ 
ing.  Firstly,  like  conventional  KS-DFT,  it  is  a  theory  that  is  exact  in  principle,  but 
practical  calculations  must  employ  an  approximate  form  for  the  unknown  exchange- 
correlation  functional.  Secondly,  unlike  conventional  KS-DFT  calculations,  the  em¬ 
bedding  formulation  introduces  an  NAKP  because  the  KS  orbitals  for  subsystem  A 
are  not  necessarily  orthogonal  to  those  of  subsystem  B.  Without  knowledge  of  the 
exact  functional  for  the  non-interacting  kinetic  energy,  this  creates  a  second  source  of 
approximation  in  the  e-DFT  approach.  The  significance  of  the  NAKP  is  system  de¬ 
pendent,  with  the  most  severe  cases  including  those  for  which  the  subsystem  densities 
Pa  and  pe  greatly  overlap. 4,9,1  ,15 

The  non-interacting  kinetic  energy  for  the  density  corresponding  to  a  set  of  A f 
closed-shell  orbitals  is 

Af 

T.{P)  =  2£<*|--V2|*).  (1.12) 

i= 1 

Standard  approximations  to  this  kinetic  energy  functional  include  the  Thomas-Fermi 
(TF)  result  for  the  homogenous  electron  gas,16,17 


Ttf  [p]  —  Ctf  /  p5//3(r)dr 


(1.13) 


where  C*tf  =  —  (37r2)2,/3,  and  the  von  Weizsacker  (vW)  result  for  the  limit  of  a 
one-electron  density, 


(1,14) 


Other  approximate  kinetic  energy  functionals  can  be  constructed  using  the  strate¬ 
gies  from  the  development  of  exchange-correlation  functionals.  For  example,  the 
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PW91k  kinetic  energy  functional1011  employs  the  analytical  form  of  the  Perdew- 
Wang  (PW91)  exchange  functional, 19  and  the  TW02  functional13  and  the  PBE2, 
PBE3,  and  PBE4  functionals9  utilizes  the  form  suggested  by  Becke.20  These  func¬ 
tionals  have  been  shown  to  successfully  describe  weakly  interacting  systems  and  co¬ 
ordination  compounds.9  Furthermore,  kinetic  energy  functionals  that  have  been  de¬ 
veloped  using  linear  response  corrections  to  the  homogeneous  electron  gas  and  have 
been  shown  to  work  well  for  metals. 3,21-23  However,  no  approximate  kinetic  energy 
functional  has  been  demonstrated  to  yield  accurate  results  for  embedded  subsystems 
that  are  connected  by  a  covalent  bond. 9,15,2 

1.3  The  Exact  Non- Additive  Kinetic  Potential 

For  each  iteration  of  the  KSCED  equations  (Eqs.  1.1  and  1.2),  {^f }  and  {0B}  (and 
thus  pa  and  pe)  are  known  from  either  the  previous  iteration  or  the  initial  guess,  and 
the  NAKP  must  be  calculated.  We  employ  a  two-step  protocol  to  obtain  the  exact 
NAKP.  In  the  first  step,  a  Levy  constrained  search25  (LCS)  or  equivalent  method  is 
used  to  determine  the  full  set  of  orthogonal  KS  orbitals,  that  correspond  to 

the  total  density  pab-  In  the  second  step,  the  NAKP  is  calculated  from  the  orbital 
sets  {</>fB},  and  {0B}. 

1.3.1  Step  1:  The  Levy  Constrained  Search 

Given  a  total  electron  density  pab,  the  fully  orthogonal  KS  orbitals  can  be 

calculated  from  an  LCS,  in  which  the  non-interacting  kinetic  energy  is  minimized  with 
respect  to  one-electron  orbitals  that  are  constrained  to  yield  Pab-25  Alternatively,  we 
employ  the  approach  of  Zhao,  Morrison,  and  Parr  (ZMP), 26-28 


in  which  the  full  set 


11 


of  KS  orbitals  are  obtained  by  solving  the  one-electron  equations 


IN nuc  ry 


v2-  y  '  -  \- 

>  ^Ir-ILI  c 


/AB,A/  x  _  /AB,A/ 


i  =  1,..., 


Nab 

~1T’ 

(1.15) 


where  1YAb  =  NA  +  NB, 


KA(r)  =  A 


PAB(r')  -  PAB(r') 


r  —  r 


dr', 


(1.16) 


Aab/2 


Pab  (r)  =  2  E  l^fB(r)|2)  and  V(,A(r)  is  a  potential  energy  function  that  restrains 
1=1 

the  Pab(i~)  to  the  target  density  Pab(i")-  Solution  of  Eq.  1.15  in  the  limit  A  — *  oo  is 
equivalent  to  performing  the  LCS. 26-28 

In  practice,  Eq.  1.15  is  solved  for  six  large,  but  finite,  values  of  A,  and  the  KS 

orbitals  and  eigenvalues  are  obtained  via  extrapolation. 26-28  For  each  value  of  A,  the 

{eA},  {0fB’A},  and  {V20fB,A}  are  calculated  and  stored  on  a  spatial  grid.  For  the 

orbitals,  extrapolation  to  A  — »  oo  is  performed  via  expansion  to  third  order  in  — , 

A 


i  AB,A/ 


r  = 


-AB 

i 


r)  +  la.a)(r)  +  ^af’(r)  +  ^5ai3,(r) 


(1.17) 


with  a  linear  least-squares  fit  of  the  expansion  coefficients  (0fB(r),  a^(r),  af\ r),  a)3^(r)} 
at  each  value  of  r.  The  { V20fB}  are  similarly  obtained  via  extrapolation  at  each  value 
of  r,  while  each  et  requires  only  a  single  extrapolation.  With  the  {4>fB}  and  {V20fB} 
obtained  on  the  spatial  grid,  the  non-interacting  kinetic  energy  for  the  total  system 
can  be  calculated  via  numerical  integration  using  Eq.  1.12. 


1.3.2  Step  2:  Exact  Kinetic  Potentials  from  KS  Orbitals 

To  calculate  the  NAKP  from  the  orbital  sets  {0fB},  {4>i~ },  and  {4>f},  we  extend  the 
approach  developed  by  King  and  Handy. 29 
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Minimization  of  the  electronic  energy  with  respect  to  the  total  electron  density 
Pab  yields  the  stationary  condition8 


STs[p] 

5p 


(r)  +  vne(r)  +  Vj[pab]  r]  +  vxc[pAB]  r] 

P=PAB 


Pab 


(1.18) 


where  pAB  is  a  Lagrange  multiplier  that  imposes  the  constraint 
Furthermore,  rearrangement  of  the  usual  KS  equations  yields 


J  Pab (r) dr  =  NAB. 


Nab/2  ,  . 

- T~\  (  -F^AB(r)V2<^B(r)  “  ei^B(r)2  )  +^e(r)  +vj[pAB]  r]  +  uxc[pab;  r]  =  0. 

PABfrj  ^  \  Z  J 

(1.19) 

Comparison  of  these  two  results  leads  to  an  exact  expression  for  the  total  kinetic 
potential, 29 


Analogous  results  can  be  derived  for  each  of  the  embedded  subsystems, 
cally,  the  electron  density  for  subsystem  A  also  obeys  a  stationary  condition, 


(1.20) 

Specih- 
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STM 


Sp 


P=PA 


(r)  +  une(r)  +  Uj[pab;  r]  +  vxc[pAB]  r]  +  unad[pA,PB;  r]  =  pA,  (1.21) 


where  pA  is  the  Lagrange  multiplier  that  imposes  the  constraint  pA(r)dr  =  Na. 
Combination  of  Eq.  1.21  with  Eq.  1.1  results  in  an  exact  expression  for  the  subsystem 
kinetic  potential, 


STs[p] 


5p 


2  1 


p=pa  PaM  “  V  2 


(r)V20f  (r)  -  4$ (r)2  )  +  pA,  (1-22) 


which  can  be  compared  with  kinetic  potential  for  the  total  system  in  Eq.  1.20. 


Insertion  of  Eqs.  1.18  and  1.21  into  Eq.  1.9  yields  pAB  =  pA ,  and  since  A  is  an 
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arbitrarily  chosen  subsystem,  we  likewise  obtain  // ab  =  hB,  or 

M  A  =  h  B-  (1-23) 


This  result  has  a  simple  physical  interpretation.  In  the  zero  temperature  limit,  the 
Lagrange  multipliers  Ha  and  /i b,  are  equivalent  to  the  chemical  potential  for  the  sub¬ 
system  electronic  densities.8  Solution  to  the  KSCED  equations  thus  yields  densities 
that  are  in  equilibrium  with  respect  to  the  number  of  electrons  in  each  subsystem. 

Finally,  insertion  of  Eq.  1.20  and  1.22  into  Eq.  1.9  yields  the  desired  expression 
for  the  NAKP, 


Note  that  the  ZMP  protocol  generally  yields  a  constant  shift  in  the  calculated  set  of 
KS  eigenenergies,  {e^};26  in  Eq.  1.24,  we  see  that  this  leads  only  to  a  constant  shift  in 
^nad  [pa >  Pbi  i"]  and  causes  no  change  in  any  calculated  observables.  Throughout  this 
study,  the  NAKP  is  shifted  such  that  it  approaches  zero  at  large  distances. 

Visscher  and  coworkers  previously  observed  that  the  NAKP  can  be  expressed  in 
terms  of  the  stationary  condition  for  the  total  system  (Eq.  1.18)  and  a  subsystem 
(Eq.  1.21). 30  They  suggested  a  strategy  in  which  the  ZMP  method  is  used  to  ob¬ 
tain  the  external  potentials  for  both  the  full  system  and  the  subsystem,  so  that  the 
difference  between  those  potentials  is  equal  to  the  NAKP,  but  this  strategy  has  not 
been  implemented  to  our  knowledge.  The  approach  presented  here  directly  expresses 
the  NAKP  in  terms  of  the  KS  orbitals  for  the  total  system  and  the  subsystem  (Eq. 
1.24);  it  avoids  performing  the  ZMP  method  for  the  subsystem,  and  it  avoids  using 
the  difference  between  two  potentials  that  are  obtained  via  the  ZMP  method.  It 
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is  straightforward  to  show  that  Eq.  1.24  recovers  the  limit  for  weakly  overlapping 
subsystem  densities  that  is  reported  in  Ref.  30. 

In  another  approach  that  does  not  utilize  the  exact  framework  of  the  KSCED 
equations,  Aguado  and  coworkers  employ  an  embedding  strategy  in  which  the  ZMP 
formalism  is  used  to  constrain  the  sum  of  the  subsystem  densities  to  that  of  the  total 
density. 31,32  This  approach  has  been  pursued  as  a  useful,  but  approximate,  strategy 
for  partitioning  a  total  density  into  subsystems. 

Other  e-DFT  strategies  also  express  the  kinetic  potential  in  terms  of  the  KS  or¬ 
bitals,  as  we  have  done  here.  For  example,  Huang  and  Carter  report  an  explicit 
expression  for  the  kinetic  potential  in  terms  of  the  KS  orbitals,  using  the  assumption 
that  the  non-interacting  kinetic  energy  is  a  linear  functional  of  the  density;  an  em¬ 
pirical  parameter  is  included  in  their  result  to  account  for  non-linear  effects.33  The 
approach  presented  here  involves  no  adjustable  parameters  and  no  assumptions  about 
the  linearity  of  the  kinetic  energy  functional. 

1.3.3  Computational  Details 

Calculations  are  performed  on  four  atomic  systems:  Li,  Ne7+,  Qjj"3"5  and  Be,  where 
Q25'5  is  a  model  3-electron  atom  that  has  a  nuclear  charge  of  +2.5.  In  all  e-DFT 
calculations,  we  take  pa  to  be  the  density  for  a  single  2s  electron,  and  pe  includes 
all  other  electrons.  The  KSCED  equations  for  each  system  were  solved  with  pe  fixed 
at  the  density  obtained  from  the  corresponding  orbitals  of  an  unrestricted  KS-DFT 
calculation  on  the  full  system;  this  is  justified  for  the  cases  studied  here  because 
solution  of  the  KSCED  equations  for  p\  subject  to  a  fixed  pe  <  po  (at  all  r),  where 
p0  is  the  exact  ground  state  density  for  the  full  system,  ensures  the  exact  calculation 
of  the  ground  state  energy  and  ground  state  density. 4  All  calculations  were  performed 
using  in-house  codes,  and  all  results  are  reported  in  atomic  units. 
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1.3. 3.1  Basis  Sets 

All  calculations  were  performed  using  the  fully  uncontracted  cc-pVTZ  basis  set  of 
Gaussian-type  orbitals  (GTOs),34  with  only  the  s-type  orbitals  included.  For  calcu¬ 
lations  on  Q^5'5,  the  Li  basis  set  was  used.  For  Ne7+,  the  most  diffuse  s-orbital  was 
removed  to  facilitate  convergence.  Although  not  reported,  all  calculations  were  also 
repeated  with  slater-type  orbitals,  which  led  to  somewhat  improved  convergence  but 
very  similar  numerical  accuracy. 

1.3. 3. 2  DFT  Implementation  Details 


Figure  1.1:  The  difference  between  the  non-interacting  kinetic  energy  Ts[p]  from 
KS-DFT  and  from  the  ZMP  method,  plotted  as  a  function  of  7.  The  extrapolation  is 
performed  using  {A}  =  {7  —  jr},  j  =  5, 4, . . . ,  0,  and  using  r  of  10  (red),  20  (green), 
and  40  (blue).  See  text  for  details. 

For  all  applications  considered  here,  p\  is  an  open  shell  system,  and  the  calcula¬ 
tions  were  performed  using  the  unrestricted  KS  formalism.  Details  for  the  unrestricted 
KSCED  equations  are  given  in  the  Appendix.  All  calculations  employ  the  Slater  ex¬ 
change  functional35  and  the  Vosko,  Wilk,  and  Nusair  correlation  functional. 36  In 
calculating  V^|SCED[pa,  Pb;  r],  a  uniform  radial  grid  is  used  to  evaluate  the  exchange- 
correlation  functions,  {0AB},  {V20AB}  ,  and  the  NAKP.  Upon  convergence  of  the 
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KSCED  equations,  the  same  radial  grid  is  used  to  evaluate  the  exchange-correlation 
energy  and  to  numerically  integrate  the  kinetic  energy.  For  Be  and  Li  the  grid  ex¬ 
tends  from  r  =  0  to  15,  while  for  Q25'5,  r  =  0  to  20  and  Ne7+,  r  =  0  to  2.  For  Be, 
Li  and  Q  the  grid  density  is  2000  points/a0  and  for  Ne7+,  20000  points/a0.  We  note 
that  future  applications  that  employ  either  a  non-uniform37  or  variational38  mesh 
will  require  fewer  gridpoints  to  achieve  the  same  level  of  accuracy.  LInless  otherwise 
stated,  the  iterative  solution  of  the  KSCED  equations  were  deemed  converged  when 
the  total  energy  of  the  system  changed  by  less  than  10~8  Hartrees  between  successive 
iterations. 

1.3. 3. 3  ZMP  Extrapolation 

To  examine  the  extrapolation  error  associated  with  the  ZMP  method,  convergence 
tests  were  performed  for  the  Li  atom  system.  The  total  density  for  the  system,  /?ab, 
and  the  reference  value  for  the  non-interacting  kinetic  energy  were  calculated  from  a 
full  KS  calculation.  This  pab  was  used  to  define  the  restraint  potential  (Eq.  1.16), 
and  the  ZMP  extrapolation  was  performed  using  six  equally  spaced  values  of  A  (i.e. , 
{A}  =  {7  —  jr},  where  j  =  5,4, .. .  ,0).  For  a  given  pair  of  parameters  7  and  r, 
the  non-interacting  kinetic  energy  was  numerically  integrated,  and  the  extrapolation 
error  was  taken  to  be  the  difference  between  this  result  and  the  reference  value  from 
the  full  KS  calculation.  Fig.  1.1  presents  this  calculated  error  as  a  function  7  and 
for  various  values  of  r.  These  results  indicate  that  the  extrapolation  error  decreases 
to  within  0.1  mH  for  7  >  500,  and  the  spacing  parameter  r  has  only  a  small  effect. 
The  error  decreases  to  within  1  pH  for  larger  values  of  7.  Results  reported  hereafter 
employ  7  =  600  and  r  =  10.  The  orbitals  from  Eq.  1.15  were  deemed  converged  when 
all  occupied  orbital  coefficients  changed  less  than  10-7  between  successive  iterations. 

The  ZMP  extrapolation  scheme  used  here  does  not  constrain  the  normalization  of 
the  orbitals.  In  general,  we  found  that  extrapolation  violated  normalization  by  less 
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then  0.01%,  and  it  was  found  that  normalizing  the  orbitals  after  extrapolation  led  to 
less  than  0.1  mff  change  in  the  total  energy.  The  results  reported  here  do  not  include 
a  posteriori  orbital  normalization. 


1.4  Results 


Figure  1.2:  The  2s  electron  density  (p\)  for  (A)  the  Q^s"5  ion,  (B)  the  Li  atom, 
(C)  the  Ne+'  ion,  and  (D)  the  Be  atom.  Calculations  performed  using  e-DFT  with 
the  non-additive  kinetic  energy  calculated  using  our  exact  protocol  (red),  the  TF 
functional  (blue),  and  the  TFvW  functional  (green).  The  black  curve,  which  is  nearly 
indistinguishable  from  the  exact  protocol,  presents  the  results  from  the  full  KS-DFT 
calculation. 


e-DFT  was  performed  for  a  series  of  three-electron  systems,  Q25'5,  Li,  and  Ne7+, 
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as  well  as  the  four-electron  Be  atom.  For  each  application,  px  was  chosen  to  include 
a  single  2s  electron,  and  the  remaining  electrons  were  included  in  pB.  In  addition 
to  using  the  exact  embedding  protocol  described  here,  the  NAKP  in  the  embedding 
calculations  was  treated  using  the  approximate  TF  kinetic  energy  functional  (Ttf  [p]  , 
Eq.  1.13)  and  the  TFvW  functional  with  the  standard  1/9  mixing  parameter  (Xtf[p]  + 
-TvW[p]). 

Fig.  1.2  presents  the  pa  obtained  in  these  e-DFT  calculations.  For  reference, 
Fig.  1.2  also  includes  the  2s  orbital  density  from  the  full  KS-DFT  calculation.  Ab¬ 
solute  agreement  between  the  KS-DFT  results  and  the  e-DFT  results  would  only 
be  expected  if  all  results  were  obtained  with  the  exact  exchange-correlation  func¬ 
tional.  Nonetheless,  since  all  calculations  in  this  study  employ  the  same  approximate 
exchange-correlation  functional,  comparison  of  the  e-DFT  and  KS-DFT  results  pro¬ 
vides  a  significant  test  of  the  accuracy  of  the  various  NAKP  descriptions. 

Fig.  1.2  clearly  demonstrates  the  sensitivity  of  e-DFT  calculations  to  the  method 
of  treating  the  NAKP.  In  comparison  to  KS-DFT,  the  e-DFT  results  from  the  approx¬ 
imate  TF  and  TFvW  functionals  are  peaked  at  significantly  shorter  radial  distances, 
and  they  qualitatively  fail  to  capture  the  nodal  structure.  Interestingly,  the  vW 
correction  to  the  TF  functional  actually  worsens  the  agreement  with  the  KS-DFT 
reference.  The  exact  embedding  protocol  describe  here,  however,  is  graphically  indis¬ 
tinguishable  from  the  KS-DFT  result. 

Further  evaluation  of  the  e-DFT  methods  can  be  obtained  by  comparing  the  cal¬ 
culated  one-electron  ionization  energies  (IEs)  for  the  various  methods.  The  e-DFT 
IE  is  calculated  from  the  difference  between  the  total  electron  energy  from  Eq.  3.6 
and  the  energy  from  a  full  KS-DFT  calculation  performed  on  the  ionized  (N-l  elec¬ 
tron)  system.  These  results  are  presented  in  Table  1.1,  which  again  illustrates  the 
qualitative  shortcomings  of  the  approximate  NAKP  treatments.  For  the  approximate 
NAKP  descriptions,  the  relative  error  between  the  e-DFT  result  and  the  KS-DFT 
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result  for  the  IEs  ranges  from  30-60%  for  3-electron  systems,  and  up  to  80%  for  Be. 
As  has  been  observed  previously, 40  including  the  vW  gradient  correction  decreases 
the  accuracy  of  the  IE  calculation.  The  exact  embedding  protocol  almost  completely 
eliminates  these  differences  with  the  reference  calculation,  with  errors  of  less  than 
0.2%  for  Q2  s'5,  Li,  and  Be  and  with  an  error  of  less  4%  for  Ne7+. 

The  lower  accuracy  of  our  embedding  protocol  for  the  case  of  Ne7+  arises  from 
the  description  of  the  nuclear  cusp.  The  KSCED  equations  converged  slowly  for  this 
case,  and  the  convergence  threshold  had  to  be  raised  to  10-5  hartrees.  By  changing 
from  GTOs  to  Slater-type  orbitals  (results  not  shown),  the  convergence  problem  was 
removed,  and  it  was  found  that  for  all  four  applications,  the  IEs  obtained  using  our 
e-DFT  protocol  were  within  1%  of  the  full  KS-DFT  result.  Below,  we  describe  how 
the  use  of  a  simple  switching  function  for  the  NAKP  in  the  cusp  region  also  removes 
these  convergence  problems  for  the  GTOs,  while  preserving  the  accuracy  of  the  IE 
calculation. 

We  note  that  the  ionization  of  the  closed  shell  Be  atom  presents  an  electronic 
structure  challenge  that  is  similar  to  the  homolytic  cleavage  of  a  covalent  bond.  From 
the  perspective  of  the  NAKP,  this  atomic  system  is  especially  challenging  since  both 
electrons  in  the  2s  “bond”  are  co-localized  on  a  single  attractive  center.  The  difficulty 
of  this  particular  case  is  confirmed  by  the  especially  poor  description  provided  by  the 
TF  and  TFvW  functionals  for  the  IE  of  the  Be  atom.  The  excellent  accuracy  of  the 
new  embedding  protocol  for  this  case  suggests  that  the  method  will  allow  for  accurate 
e-DFT  calculations  in  which  the  subsystems  are  linked  by  covalent  bonds. 

Fig.  1.3  illustrates  the  KSCED  potentials,  KefpCED[pA,  Pb;  r],  and  the  correspond¬ 
ing  NAKPs,  tw[dA,PB;  r],  that  are  obtained  from  the  exact  embedding  calculations. 
For  each  system,  the  similarity  between  these  two  potentials  illustrates  the  domi¬ 
nance  of  the  NAKP  at  short  distances.  However,  the  NAKP  decays  rapidly,  and  the 
KSCED  potential  is  dominated  at  larger  distances  by  the  Coulombic  terms  (Eq.  1.5). 
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Table  1.1:  Total  energy  (TE)  and  ionization  energy  (IE)  obtained  using  KS-DFT 
and  e-DFT. 


Atom 

Calculation 

TE 

IE 

Error  (%) 

Qzf 

KS 

-4.799405 

0.060142 

TFvW 

-4.835122 

0.095258 

59.99 

TF 

-4.819221 

0.079357 

33.28 

Exact  Embedding 

-4.799510 

0.060247 

0.18 

Li 

KS 

-7.343870 

0.201098 

TFvW 

-7.443321 

0.300549 

49.45 

TF 

-7.408509 

0.265737 

32.14 

a  Exact  Embedding 

-7.344046 

0.201274 

0.09 

Ne7+ 

KS 

-101.964612 

8.754056 

TFvW 

-106.413890 

13.203334 

50.83 

TF 

-105.630042 

12.419486 

41.87 

b  Exact  Embedding 

-102.294207 

9.083650 

3.77 

Be 

KS 

-14.447017 

0.331698 

TFvW 

-14.717243 

0.601924 

81.47 

TF 

-14.635950 

0.520631 

56.96 

a  Exact  Embedding 

-14.447463 

0.328900 

0.13 

a  KSCED  equations  converged  to  10  '  Hartree 
b  KSCED  equations  converged  to  10-5  Hartree 
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Although  it  is  not  visible  from  the  scale  of  the  plots  in  Fig.  1.3,  the  unad[PA,  Pb;  r] 
term  comprises  less  than  1%  of  the  f/gSCED[pA,  Pb;  r]  for  distances  greater  than  3  a.u. 
for  all  cases.  (For  Ne7+,  this  regime  is  reached  at  0.43  a.u.) 
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Figure  1.  3:  The  KSCED  effective  potential,  14esced[pa,  Pb;  r],  for  (A)  the  Q2g5 
ion,  (B)  the  Li  atom,  (C)  the  Ne+'  ion,  and  (D)  the  Be  atom  and  the  NAKP, 
hmd[pA,  Pb;  r],  for  (E)  the  Q^°'5  ion,  (F)  the  Li  atom,  (G)  the  Ne+'  ion,  and  (H) 
the  Be  atom  using  the  e-DFT  protocol  presented  here. 


Comparison  of  the  NAKPs  in  Fig.  1.3E-H  with  the  densities  in  Fig.  1.2  illustrates 
that  the  nodal  structure  in  the  2s  electron  density  is  enforced  by  the  NAKP.  For  each 
system,  the  large  outer  peak  in  the  NAKP  coincides  with  the  nodal  feature  in  the 
2s  density.  Unlike  the  KS-DFT  results,  we  note  that  the  densities  obtained  using 
e-DFT  in  Fig.  1.2  do  not  exhibit  a  genuine  radial  node,  since  p\  corresponds  to  the 
ground-state  eigenvector  of  Eq.  1.1.  A  large  peak  in  the  e-DFT  effective  potential 
is  therefore  essential  to  achieve  the  correct  2s  shell  structure.  The  NAKPs  obtained 
from  the  approximate  TF  and  TFvW  functionals  do  not  exhibit  this  pronounced  peak 
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(not  shown),  which  leads  to  the  poor  descriptions  for  the  2s  electron  density  (Fig. 
1.2)  and  the  IE  (Table  1.1). 

In  addition  to  the  pronounced  outer- most  peak  for  each  NAKP  in  Fig.  1.3E-H, 
oscillations  at  short  distances  are  observed.  This  oscillatory  behavior  is  sensitive  to 
the  basis  set  representation.  Small  changes  in  the  orbital  coefficients  for  regions  of  low 
density  give  rise  to  large  changes  in  the  kinetic  potential  (Eq.  1.24),  resulting  in  slow 
convergence  of  the  KSCED  equations.  (These  oscillations  are  not  observed  when  the 
density  vanishes  at  large  distances  since  the  basis  set  expansion  is  dominated  by  only 
the  slowest-decaying  function  in  that  regime.)  Using  STOs  rather  than  GTOs,  the 
NAKP  oscillations  at  short  distances  were  diminished  (not  shown),  and  the  iterative 
convergence  was  improved.  In  future  applications  of  the  exact  embedding  protocol 
with  GTOs,  the  use  of  the  convergence  acceleration  algorithms  such  as  DIIS39  may 
prove  beneficial.  However,  we  now  demonstrate  that  the  problems  associated  with 
poor  convergence  and  NAKP  oscillations  can  be  alleviated  with  a  simple  modification 
of  Eq.  1.24. 

As  pa  vanishes  close  to  the  nucleus,  evaluation  of  the  second  term  in  Eq.  1.24  be¬ 
comes  unstable,  leading  to  the  aforementioned  convergence  problems.  This  is  avoided 
by  introducing  a  switching  function  that  changes  from  the  exact  expression  for  the 
kinetic  potential  of  subsystem  A  to  the  corresponding  TF  approximation  near  the 
nucleus, 

o  Nab/2  f  i  \ 

V„ad[pA,PB;r]  =  - -  (L25) 

PAB(r)  “  V  Z  J 

Na/ 2  , 

-  yyy  E  -  y^w2)  x  (1.26) 

(1  - /[pb;  <-])  -  Qctf  /'A  3)  / [pb ;  r] , 
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where  f[pv]r]  is  the  smooth  switching  function 

/[PB;r]  =  eK(_pB(r)+^}  +  1-  (1.27) 

Previous  work  used  a  similar  function  to  switch  between  approximate  expressions  for 
the  NAKP  in  the  vicinity  of  the  nuclear  cusp.40  The  parameters  p'B  and  k,  determine 
the  radial  distance  and  the  abruptness  with  which  switching  occurs,  respectively.  The 
parameter  p'B  was  related  to  the  integrated  electron  density  in  the  cusp  region,  setting 
Pb  =  Pb  (r/),  where 

prf 

£  =  47T  /  r2pB(r)dr.  (1-28) 

Jo 

e-DFT  results  obtained  using  range  of  values  for  k  and  £  were  compared  to  deter¬ 
mine  robust  parameters  for  the  switching  function.  Setting  n  =  50,  we  varied  £  over 
the  range  from  0.4  to  0.8  for  Li  and  Ne7+,  which  led  to  changes  in  the  total  calculated 
energy  of  less  than  0.4  mH  and  5  mH,  respectively.  Similarly,  setting  £  =  0.6  and 
varying  k  over  the  range  from  50  to  500  led  to  differences  of  less  than  0.1  mH  for 
both  Li  and  Ne. 

Using  the  NAKP  expression  in  Eq.  1.27  with  £  =  0.6  and  k  =  50,  our  e-DFT  pro¬ 
tocol  was  applied  to  all  four  systems,  and  the  results  are  presented  in  Table  1.2.  All 
calculations  reached  full  10~8  convergence  within  80  iterations  of  the  KSCED  equa¬ 
tions,  in  contrast  with  the  calculations  using  Eq.  1.24,  which  was  difficult  to  converge 
in  some  cases  even  with  2000  iterations.  Furthermore,  the  e-DFT  calculations  with 
the  modified  NAKP  expression  in  Eq.  1.27  yields  excellent  accuracy  in  comparison 
to  the  full  KS-DFT  equations,  with  less  than  1.5%  error  in  the  IE  for  all  cases. 

For  the  Li  atom,  Fig.  1.4  compares  the  NAKP,  the  KSCED  effective  potential, 
and  the  2s  electron  density  obtained  by  solving  the  KSCED  equations  using  either 
Eq.  1.24  (black)  or  Eq.  1.27  (red)  for  the  NAKP.  The  black  curves  in  this  figure  are 
the  same  as  those  for  Li  in  Figs.  1.2  and  1.3.  It  is  clear  from  Fig.  1.4A  that  at  short 
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Table  1.2:  Total  energy  (TE)  and  ionization  energy  (IE)  obtained  using  e-DFT  with 
the  NAKP  switching  function  (Eq.  1.27). 


Atom 

TE 

IE 

Error  (%) 

(y-°-5 

M22.5 

-4.799142 

0.059879 

0.44 

Li 

-7.342720 

0.199948 

0.57 

Ne7+ 

-101.843497 

8.632941 

1.38 

Be 

-14.443703 

-0.328383 

1.00 

distances,  the  switching  function  produces  a  relatively  featureless,  repulsive  NAKP 
due  to  the  TF  approximation;  the  arrow  in  this  figure  indicates  the  radial  distance  r' 
that  corresponds  to  the  parameter  £  =  0.6.  Fig.  1.4B  illustrates  that  the  repulsive 
NAKP  largely  cancels  the  attractive  electron- nuclear  Coulomb  term  in  the  KSCED 
effective  potential  (Eq.  1.5).  As  p\  vanishes  at  the  nucleus,  the  KSCED  effective 
potential  must  also  approach  zero.30  The  remaining  oscillations  at  radial  distances 
in  Fig.  1.4B  are  an  artifact  of  the  finite  basis  set.  Finally,  Fig.  1.4C  demonstrates 
that  the  2s  electron  density  that  is  obtained  using  the  switching  function  does  not 
reproduce  the  features  of  the  radial  node,  but  it  recovers  the  exact  embedding  result 
for  distances  beyond  1  a.n.  This  close  agreement  at  large  distances  is  expected41 
from  the  accuracy  of  the  IE  calculations  in  Table  1.2.  In  light  of  the  much  improved 
convergence  efficiency,  use  of  the  NAKP  expression  in  Eq.  1.27  compares  favorably 
with  exact  embedding  via  Eq.  1.24. 

1.5  Extension  to  Larger  Systems 

The  calculations  reported  here  demonstrate  a  proof-of-principle  for  the  exact  calcula¬ 
tion  of  the  NAKP.  However,  for  applications  of  e-DFT  to  large  systems,  performance 
of  the  ZMP  extrapolation  at  each  iteration  of  the  KSCED  equations  is  impractical. 
Nonetheless,  the  short-ranged  nature  of  the  NAKP  (see  Fig.  1.3E-H)  suggests  several 
strategies  for  employing  our  e-DFT  protocol  in  larger  systems. 

For  example,  suppose  that  subsystem  B  is  further  divided  into  fragments  (Bl5  B2, 
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Figure  1.4:  (A)  The  NAKP,  (B)  the  KSCED  effective  potential,  and  (C)  the  2s 
electron  density  (/?a)  for  the  Li  atom,  obtained  using  exact  embedding  (black)  and 
using  the  modified  NAKP  in  Eq.  1.27  (red).  The  arrow  indicates  the  radial  distance 
at  which  switching  occurs. 


. . .,  Bj),  and  consider  the  sum  of  the  NAKP  terms  due  to  the  individual  fragments, 


^nad  [PA  i  Pli  i  f] 


/ 

£ 

i—  1 


(sr.\p\ 

ST.\p] 

\  *P 

P=PA+PBi  Sp  P 

P=PA 


(1.29) 


This  equation  is  exact  in  the  limit  of  one  fragment,  and  its  implementation  with  our 
protocol  will  avoid  ZMP  extrapolation  for  anything  larger  than  the  union  of  subsystem 
A  with  a  single  fragment. 

The  assumption  in  Eq.  1.29  that  the  NAKP  is  additive  among  the  fragments 
must  be  tested.  However,  any  error  introduced  from  this  assumption  can  be  partially 
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corrected  using  the  TF  (or  other)  approximate  kinetic  energy  functional, 


Wad  [PA  i  PB  •  I". 
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(1.30) 


Here,  the  first  term  on  the  RHS  corresponds  to  the  NAKP  obtained  from  the  TF 
functional  for  the  full  system.  In  the  second  term,  the  contribution  due  to  each  of  the 
fragments  using  the  TF  approximation  is  removed,  and  in  the  third  term,  each  of  the 
fragment  contributions  is  replaced  using  the  exact  protocol.  The  short-ranged  nature 
of  the  NAKP  suggests  that  distance-based  cutoffs  can  be  employed  with  summations 
in  Eqs.  1.29  and  1.30,  allowing  for  significant  computational  savings. 


1.6  Conclusions 

We  have  described  a  general  and  exact  protocol  for  treating  the  non-additive  kinetic 
potential  in  embedded  density  functional  theory  calculations.  In  applications  to  a 
series  of  three-  and  four-electron  systems,  we  have  numerically  demonstrated  the  ap¬ 
proach,  and  we  have  illustrated  the  qualitative  failures  that  can  arise  from  the  use 
of  approximate  kinetic  energy  functionals.  We  have  also  shown  that  improved  con¬ 
vergence  of  the  KSCED  equations  can  be  obtained  with  appropriate  switching  of  the 
NAKP  in  the  vicinity  of  the  nuclear  cusps,  and  we  have  described  possible  strategies 
for  the  scalable  implementation  of  our  embedding  protocol  in  large  systems.  Natu¬ 
ral  applications  of  exact  embedding  include  the  rigorous  calculation  of  one-electron 
pseudopotentials,  the  calculation  of  DFT  embedding  potentials  for  use  with  high-level 
ab  initio  calculations  on  small  subsystems, 5,21,42,43  and  the  accurate  implementation 
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of  the  “molecular  embedding”  strategy  in  which  each  molecule  of  a  large  system  is 
assigned  to  a  different  embedded  subsystem. 44 


1.7  Appendix:  Unrestricted  Open-Shell  e-DFT 

For  unrestricted  open-shell  e-DFT  calculations,  the  density  of  each  subsystem  is  fur¬ 
ther  partitioned  into  a  and  /?  spin  densities,  such  that  pab  =  Pa  +  Pa  +  Pb  +  Pb-  This 
leads  to  the  KSCED  equations 
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Here,  N"  is  the  number  of  electrons  in  each  subsystem,  and  p^(r)  =  Eir'wi2, 
where  /i  €  {A,  B}  and  u  e  {a,  /?}.  The  KSCED  effective  potential,  KgSCED[p^,  p^,  p^,  p^;  r], 

is 


rJSCED[PA.  Pa,  Pb,  Pb!  rl  =  ivM+pjIpab;  i-J+ivKpa+Pb),  (pKpb);  r]-HWd[pA,  Pb;  r] 

(1.35) 

where  une(r)  and  uj[pab;  r]  are  unchanged  from  Eq.  1.5,  uxc[(pa  +  Pb),  (Pa  +  Pb);  r] 
is  the  usual  open-shell  exchange-correlation  potential  for  the  total  system,8  and  the 
NAKP  is  discussed  below. 

The  kinetic  energy  functional  is  separable  into  two  different  spin  contributions8 


r,[p“,pj]  =  r,[p“  o]  +  r,[o,p«], 


(1.36) 


where 


r.(p“  o]  =  ^(0H-|v2IC'“)  (i-37) 

i=l 

and  likewise  for  Ts[0, /cr].  Therefore,  the  NAKP  depends  only  on  spin  densities  cor¬ 
responding  to  the  same  spin,  such  that 
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The  ZMP  extrapolation  is  used  to  calculate  the  KS  spin  orbitals  and 

eigenvalues  for  the  full  system,  exactly  as  is  described  in  the  text,  except  that 

the  total  spin  density  is  employed  instead  of  the  total  electron  density.  Finally,  our 
exact  expression  for  the  NAKP  for  open-shell  systems  is  modified  from  Eq.  1.24  as 
follows: 
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The  TF  approximation  for  the  non-additive  kinetic  energy  in  an  open-shell  calcu¬ 
lation  is 


rTn? K.rfi]  =  22/3Ctf  J  (pAB/3(r)  -  ftAW  -  Pb5/SM)  dr,  (1.41) 


and  corresponding  result  for  the  TFvW  functional  is 
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Chapter  2 

Embedded  density  functional  theory  for  covalently 
bonded  and  strongly  interacting  subsystems 

2.1  Introduction 

Important  methodological  challenges  in  electronic  structure  theory  include  the  long- 
timescale  simulation  of  ab  initio  molecular  dynamics  and  the  seamless  combination  of 
high-  and  low-level  electronic  structure  methods  in  complex  systems.  Methods  that 
exploit  the  intrinsic  locality  of  molecular  interactions  have  demonstrated  encouraging 
progress  towards  these  goals. 1-17 

In  particular,  orbital-free  embedded  DFT  (e-DFT)  offers  a  formally  exact  ap¬ 
proach  to  electronic  structure  theory  in  which  the  interactions  between  subsystems 
are  evaluated  in  terms  of  their  electronic  densities. 1_  Recent  work  has  demonstrated 
that  constructing  the  embedded  subsystems  from  individual  molecules  leads  to  a 
linear-scaling  electronic  structure  approach  that  maps  naturally  onto  distributed- 
memory  parallel  computers, 13,18  and  it  provides  a  systematic  framework  for  calculat¬ 
ing  electronic  excited  states  in  condensed  phase  systems. 19,20  However,  approximate 
treatments  of  the  non-additive  kinetic  potential  (NAKP)  limit  the  accuracy  of  this 
approach  in  applications  involving  strongly  interacting  subsystems. 21  For  example, 
severe  artifacts  in  the  structure  of  liquid  water,  including  the  complete  absence  of  a 
second  peak  in  the  oxygen-oxygen  radial  distribution  function,  have  been  predicted 
from  existing  approximations  to  the  NAKP, 18  and  e-DFT  applications  involving  co- 
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valently  bonded  embedded  subsystems  have  also  been  shown  to  qualitatively  fail. 21-23 
The  development  of  improved  methods  to  address  the  NAKP  problem  will  open  new 
doors  for  on-the-fly,  massively  parallel  electronic  structure  calculations  in  general, 
condensed-phase  systems. 

In  this  paper,  we  describe  progress  towards  the  development  of  accurate,  scalable 
treatments  for  the  NAKP  in  e-DFT.  We  provide  the  first  molecular  applications  of 
our  recently  developed  Exact  Embedding  (EE)  method,24  demonstrating  that  it  suc¬ 
cessfully  describes  the  breaking  of  covalent  bonds  and  hydrogen  bonds  with  chemical 
accuracy.  Additionally,  we  introduce  and  numerically  demonstrate  a  pairwise  ap¬ 
proximation  to  the  NAKP,  which  allows  for  the  scalable  implementation  of  the  EE 
method  in  large  systems.  Benchmark  calculations  are  presented  for  systems  with  up 
to  125  molecules,  demonstrating  that  parallel  implementation  of  the  method  enables 
constant  system-size  scaling  of  the  wall-clock  calculation  time. 


2.2  Theory 


2.2.1  Orbital-Free  Embedded  DFT 

We  utilize  the  orbital- free  e-DFT  formulation  of  Cortona1  and  Wesolowski  and  cowork¬ 
ers.2’3  For  the  case  in  which  the  total  electronic  density  Pab  is  partitioned  into  two 
subsystems,  pab  =  Pa  +  Pb,  the  corresponding  one-electron  orbitals  obey  the  Kohn- 
Sham  Equations  with  Constrained  Electron  Density  (KSCED),3 


--V2  +  veff[pA,PAB;r] 
~ -V2  +  ueff[pB,  Pab!  r] 


(2.1) 

(2.2) 
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where  i  —  1, , . . ,  NA,  j  —  1, ,  1VB,  and  NA  and  NB  are  the  number  of  electrons  in 
the  respective  subsystems.  veff  is  the  effective  potential  for  the  coupled  one-electron 
equations,  such  that 


ves[pA,  Pab ;  r]  =  une(r)  +  vj[pAB]  r]  +  vxc[pAB]  r] 

"itiad  [p a  •  Pab  i  r] , 


(2.3) 


where  the  Nnuc  nuclei  occupy  positions  {R;}, 
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(2.5) 

(2.6) 


and  Exc  [p]  is  the  exchange-correlation  functional.  unad  [pa  ,  Pab  ;  r]  is  the  potential  due 
to  the  non-additive  kinetic  energy  for  non-interacting  electrons,  such  that 


^nad  [pa  j  Pab  i  J”. 


hrsnad[pA,pB] 

5pA 


(2.7) 


where  Tsnad[pA,PB]  =  Ts[pab]  —  Ts[pa]  —  Ts[pr].  The  subsystem  densities  are  con- 

na 

structed  from  the  corresponding  KS  orbitals,  using  Pa(i~)  =  Ei^mi2  and  = 

2=1 

nb 

|0B(r)|2.  Eqs.  3. 1-2.7  are  easily  generalized  for  the  e-DFT  description  of  multiple 
i= i 

embedded  subsystems. 1,18 

Two  aspects  of  e-DFT  are  worth  emphasizing.  Firstly,  like  conventional  Kohn- 
Sham  (KS)-DFT,  it  is  a  theory  that  is  exact  in  principle,  but  practical  calculations 
must  employ  approximations  to  the  unknown  exchange-correlation  functional.  Sec¬ 
ondly,  unlike  conventional  KS-DFT  calculations,  the  embedding  formulation  intro- 
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duces  the  NAKP,  unad[PA,  Pab;  r],  since  the  one-electron  orbitals  for  subsystem  A 
are  not  necessarily  orthogonal  to  those  of  subsystem  B.  Without  knowledge  of  the 
exact  functional  for  the  non-interacting  kinetic  energy,  this  creates  a  second  source 
of  approximation  in  e-DFT  calculations.  The  significance  of  the  NAKP  is  system- 
dependent,  with  the  most  severe  cases  including  those  for  which  the  subsystem  den¬ 
sities  greatly  overlap;  no  approximate  kinetic  energy  functional  has  been  previously 
demonstrated  to  yield  accurate  results  for  embedded  subsystems  that  are  connected 
by  covalent  bonds. 3’21,22,25’26 


2.2.2  Exact  Calculations  of  NAKP 

We  have  recently  developed  the  Exact  Embedding  (EE)  method  to  calculate  the 
NAKP. 21  The  general  method  can  be  summarized  for  two  embedded  subsystems  as 
follows:  A  Levy  constrained  search  (LCS)2'  or  equivalent  technique  is  first  used  to 
determine  the  full  set  of  orthogonal  KS  orbitals,  {0fB},  that  correspond  to  the  total 
density  /?ab  from  the  latest  iteration  of  Eqs.  3. 1-2. 3.  Then,  from  the  KS  orbitals 
{(/•f'5},  {0f },  and  {4>f },  the  corresponding  kinetic  potentials  are  calculated  using  the 
exact  result  of  King  and  Handy, 28 


vts  (r) 


EjLi(-%Mr)v2Mr)) 

p(r) 


(2.8) 


where  n  is  the  number  of  occupied  orbitals,  e*  is  the  KS  eigenvalue  corresponding  to 
orbital  0j,  and  /x  is  a  constant.  Finally,  the  NAKP  needed  for  the  next  iteration  of 
Eqs.  3. 1-2.3  is  calculated  directly  from  the  difference 


fn«i[PA,  Pab;  r]  =  t>r,B(r)  -  t>£,(r), 


(2.9) 
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where  the  superscripts  in  this  equation  indicate  the  orbital  set  to  which  each  kinetic 
potential  corresponds. 

Rather  than  explicitly  performing  the  LCS,  we  use  the  equivalent  protocol  of  Zhao, 
Morrison,  and  Parr  (ZMP) 29-31  to  obtain  the  exact  non-interacting  kinetic  energy  and 
the  KS  orbitals  This  requires  solution  of  the  following  one-electron  equations 


-V'  +  Uext(r)  +  vc  (r) 


<^aBM  =  e^AB(i 


(2.10) 


in  the  limit  A  — »  oo,  where  i  =  1, ,  ( NA  +  NB),  and 


A 


p(r')  -pAB(r)ir, 


(2.11) 


UeXt(r)  corresponds  to  any  well-behaved  external  potential,30’31  and  various  choices 
for  this  potential  are  described  in  Sec.  Ill  B.  In  practice,  Eq.  2.10  is  solved  at  several 
large,  finite  values  of  A,  and  the  KS  orbitals  and  eigenvalues,  as  well  as  the  final  non¬ 
interacting  kinetic  energy,  are  obtained  via  extrapolation.29  31  In  Sec.  V,  we  discuss 
a  technique  to  robustly  implement  the  ZMP  step  for  NAKP  calculations  in  large 
systems. 

The  EE  method  outlined  in  Eqs.  2.8  -  2.11  is  unique  in  that  it  allows  for  the  for¬ 
mally  exact  calculation  of  the  total  electronic  density  within  the  e-DFT  framework, 
using  integer  orbital  occupancies  and  without  approximations  to  the  NAKP.  The 
method  was  previously  demonstrated  for  atomic  systems  with  strongly  overlapping 
subsystem  densities,24  and  the  current  paper  presents  its  first  molecular  applications. 
We  note  that  several  other  groups  have  also  used  density  inversion  techniques  to  cal¬ 
culate  the  NAKP,  assuming  that  the  total  electron  density  is  already  available  from 
another  electronic  structure  calculation.  ’  ’  In  particular,  Visscher  and  coworkers 
have  applied  this  approach  to  molecular  systems  with  the  aim  of  developing  improved 
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non-additive  kinetic  energy  functionals. 23  Furthermore,  Partition  DFT  has  been  in¬ 
troduced  as  a  formally  exact  embedding  scheme  in  which  subsystem  densities  are  de¬ 
scribed  using  partially  occupied  orbitals,  and  it  has  been  applied  to  one-dimensional 
model  systems.6 

2.3  Implementation  Details 

We  have  implemented  e-DFT  in  the  Molpro  quantum  chemistry  package,31  allowing 
for  calculation  of  the  NAKP  with  either  approximate  functionals  or  the  EE  method. 
In  this  section,  methodological  and  numerical  aspects  of  the  implementation  are  dis¬ 
cussed. 

2.3.1  Supermolecular  vs.  Monomolecular  Basis  Sets 

The  atom-centered  basis  sets  used  to  solve  the  KSCED  (Eqs.  3.1  and  3.2)  are  im¬ 
plemented  using  two  different  conventions.21,35  In  the  monomolecular  basis  set  con¬ 
vention,  the  density  for  each  embedded  subsystem  is  described  using  only  the  basis 
functions  that  are  centered  on  atoms  belonging  to  that  subsystem.  In  the  super- 
molecular  basis  set  convention,  the  density  for  each  embedded  subsystem  is  described 
using  the  same  basis  set,  which  includes  functions  that  are  centered  on  all  atoms  in 
the  system.  The  supermolecular  basis  set  convention  provides  a  closer  approximation 
to  the  complete  basis  set  limit,  although  it  is  more  costly. 

2.3.2  ZMP  Step 

In  our  implementation,  the  ZMP  step  of  the  EE  method  is  performed  by  solving 
Eq.  2.10  for  six  large,  finite  values  of  A.  The  KS  orbitals  {0fB}  are  then  obtained 
from  extrapolation  of  the  atomic  orbital  coefficients  for  the  {0fB},  using  a  third- 
order  polynomial  in  A-1,  and  normalization  of  the  extrapolated  orbitals  is  enforced 
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a  posteriori.  The  KS  eigenvalues  {efB}  are  similarly  obtained  from  extrapolation  of 
the  {efB}.  Ts[pab]  is  calculated  analytically  from  the  extrapolated  orbital  coefficients, 
which  ensures  that  the  total  energy  from  the  EE  method  is  bound  from  below  by  the 
KS-DFT  energy. 

In  the  limit  A  — »  oo,  the  solutions  to  Eq.  2.10  are  independent  of  the  choice 
of  external  potential  next(r), 29-31  although  vext(r)  does  effect  the  convergence  with 
increasing  A.  Various  options  where  thus  considered,  including 


Vext(r) 

=  Aie(r), 

(2.12) 

Vext(r) 

=  We(r)+(l  iyA  +  ArB)ffi[dAB;r], 

(2.13) 

^ext(r) 

=  vne  (r)  +  vj  [ pab  ;  r]  +  vxc  [pA b ;  r] . 

(2.14) 

At  every  iteration  of  the  KSCED,  these  versions  of  uext(r)  are  all  available  without 
the  need  for  additional  computation.  Test  calculations  have  indicated  that  the  exter¬ 
nal  potential  in  Eq.  2.14  leads  to  the  fastest  convergence  of  the  extrapolation  with 
increasing  A,  and  this  potential  is  used  in  all  results  for  the  EE  method  reported  here. 

2.3.3  NAKP  Numerics  for  Regions  of  Weak  Density  Overlap 

Numerical  evaluation  of  the  kinetic  potential  from  Eq.  2.8  is  unstable  in  regions 
for  which  the  corresponding  density  vanishes.  The  problem  is  exacerbated  by  the 
incorrect  distance  dependence  of  the  low-density  tails  obtained  from  calculations  using 
Gaussian-type  orbitals  (GTOs).28  However,  these  numerically  treacherous  regions 
correspond  to  weak  overlap  between  subsystem  densities,  where  the  magnitude  of  the 
NAKP  is  necessarily  small  and  easily  approximated. 2  We  thus  utilize  a  density-based 
criterion  to  switch  from  the  exact  expression  for  the  kinetic  potential  to  a  numerically 
stable  approximation,  such  as  the  Thomas- Fermi  (TF)  kinetic  potential.  The  protocol 
used  to  perform  this  switching  is  described  below. 
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In  a  first  step,  we  calculate  the  constant  shift  that  is  needed  to  match  the  exact  re¬ 
sult  for  each  kinetic  potential  to  the  corresponding  TF  result  in  a  prescribed  switching 
region.  Specifically,  for  each  of  the  kinetic  potentials  (i.e.,  vts(  r)  G  {u^B(r),uAs(r),uBs(r)} 
which  correspond,  respectively,  to  p(r)  G  {pabM, PaM, PbM}),  the  average  differ¬ 
ence  (A  G  {Aab,Aa,Ab})  between  the  results  from  Eq.  2.8  and  from  the  TF  func¬ 
tional  is  evaluated  in  the  vicinity  of  the  p(r)  =  p  density  isosurface.  Each  A  is 
computed  over  gridpoints  in  the  region  £  <  /[p;  r]  <  (1  —  £),  where 


f\p\  r 


1 

eK(p(r)-p')  _|_  l  ’ 


(2.15) 


£,  k,  and  p  are  parameters  that  define  the  switching  region,  and  the  relative  contri¬ 
bution  from  each  gridpoint  is  weighted  according  to 


u\p]  r]  =  e“K(pAB(r)_p(r)). 


(2.16) 


Note  that  the  weighting  function  in  Eq.  2.16  is  uniform  for  the  case  of  p  =  Pab,  and 
it  preferentially  selects  values  which  p(r)  ~  Pab(f)  for  the  cases  in  which  p(r)  is  one 
of  the  subsystem  densities. 

In  a  second  step,  each  kinetic  potential  is  computed  on  the  grid;  this  is  done 
by  vertically  shifting  the  exact  result  with  the  corresponding  A  and  then  smoothly 
switching  to  the  TF  result  at  densities  below  p',  using  the  density-based  switching 
function  /[p;r]  in  Eq.  2.15.  Finally,  the  NAKP  is  calculated  from  the  smoothly 
switched  kinetic  potentials  using  Eq.  2.9.  The  vertical  shifts  that  are  applied  to 
kinetic  potentials  simply  give  rise  to  an  additive  constant  in  the  final  NAKP,  which 
has  no  physical  effect.  Although  we  find  that  switching  to  the  TF  functional  at 
low  densities  is  both  convenient  and  accurate,  the  protocol  described  above  could  be 
performed  using  any  approximate  kinetic  energy  functional. 
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2.4  Results:  Small  Systems 

2.4.1  Calculation  Details 

In  this  section,  e-DFT  calculations  are  presented  for  the  dissociation  curves  of  (H20)2 
and  the  covalently  bound  Li+-Be  and  CH3-CF3  molecules;  standard  KS-DFT  calcula¬ 
tions  are  included  for  comparison.  All  results  are  obtained  using  the  Molpro  quantum 
chemistry  package,3  with  KS-DFT  available  in  the  standard  version  and  with  the 
e-DFT  method  implemented  in  our  modified  version.  In  the  e-DFT  calculations, 
the  NAKP  is  described  using  either  the  EE  method  or  the  approximate  TF36,37  and 
LC9438  kinetic  energy  functionals;  these  approaches  will  hereafter  be  referred  to  as 
e-DFT-EE,  e-DFT-TF,  and  e-DFT-LC,  respectively. 

All  calculations  in  this  section  are  performed  using  the  B88-P86  exchange-correlation 
(XC)  functional. 39,40  Both  the  XC  functional  and  the  NAKP  are  evaluated  on  a  grid  of 
Becke-Voronoi41  cells  with  resolution  to  limit  the  integration  error  of  Slater  exchange 
to  1CT12  Hartree;  the  grid  is  generated  using  the  Molpro  directive  GRID=10-12. 

The  KSCED  in  Eqs.  3. 1-3.2  are  initialized  from  the  gas  phase  density  of  each 
subsystem,  and  the  eigensolutions  for  each  set  of  equations  are  updated  at  every 
iteration.  Convergence  of  these  equations  is  improved  with  the  molecular  orbital 
(MO)  shifting  and  direct  inversion  of  iterative  subspace  (DIIS)  algorithms.42,  3  For 
the  water  dimer,  an  MO  shift  of  -0.5  Hartree  is  employed,  whereas  a  -1.0  Hartree 
shift  is  used  for  Li+-Be  and  CH3-CF3.  Since  the  DIIS  algorithm  leads  to  slow  final 
convergence,44  it  is  discontinued  once  the  root  mean  squared  difference  (RMSD)  of 
total  density  matrix  elements  changes  by  less  than  5  x  10~4  between  two  successive 
iterations.  The  KSCED  equations  are  deemed  converged  when  the  total  energy  of  the 
system  changes  by  less  than  10-6  Hartree  and  the  RMSD  in  the  total  density  matrix 
is  smaller  than  10-5  between  two  successive  iterations. 
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For  the  ZMP  step,  extrapolation  of  the  solutions  to  Eq.  2.10  is  performed  using 
A  =  7  +  rj,  where  j  =  0, 1, . . . ,  5.  Unless  otherwise  noted,  calculations  for  the 
water  dimer  and  Li+-Be  employ  7  =  5000  and  r  =  100,  whereas  calculations  for 
CH3-CF3  employ  7  =  100  and  r  =  10.  To  reach  adequate  convergence,  Eq.  2.10  is 
solved  in  several  stages.  Firstly,  a  coarse  solution  is  reached  by  using  an  MO  shift 
of  — 103  Hartree  and  a  value  of  A  =  100.  Subsequently,  using  this  coarse  solution  as 
a  starting  point,  the  Eq.  2.10  solved  using  a  smaller  MO  shift  of  —84  Hartree  and 
with  A  =  7.  Finally,  solution  of  Eq.  2.10  for  each  increasing  value  of  A  needed  for 
extrapolation  employs  the  solution  for  the  prior  value  of  A  as  a  starting  point.  The 
D1IS  algorithm  is  used  throughout.  The  orbitals  from  Eq.  2.10  are  deemed  converged 
when  the  RMSD  in  the  density  matrix  was  smaller  than  1CT9  between  two  successive 
iterations;  significantly  tighter  convergence  is  needed  for  the  ZMP  equations  than  for 
the  KSCED,  to  ensure  an  accurate  extrapolation. 

Calculations  for  the  water  dimer  variously  employ  the  aug-pc-3,  aug-pc-2,  and 
aug-pc-1  basis  sets,45  in  each  case  using  only  the  s-  and  p-type  functions  for  the 
hydrogen  atoms  and  the  s-,  p-,  and  d-type  functions  for  the  oxygen  atoms.  These 
water  dimer  basis  sets  are  hereafter  referred  to  as  the  modified  aug-pc-3,  aug-pc-2, 
and  aug-pc-1  basis  sets,  respectively.  In  all  calculations  for  Li+-Be  and  CH3-CF3, 
the  Li,  Be,  and  C  atoms  are  described  using  the  s-,  p-,  and  d-type  functions  of  the 
combined  aug-pc-4  and,  for  Li  and  Be,  the  cc-pVQZ  (core/valence),  for  C,  the  cc- 
pV6Z  (core/valence)  basis  sets,46  and  the  H  and  F  atoms  are  described  using  the  full 
aug-pc-1  basis  set.45  Sensitivity  of  the  e-DFT  calculations  to  the  basis  set  is  discussed 
in  the  next  section. 

Larger  basis  sets  provide  a  better  description  of  low-density  regions,  allowing  for 
the  use  of  smaller  values  for  the  parameter  p  in  Eqs.  2.15  and  2.16  and  providing 
robustness  with  respect  to  the  choice  of  this  parameter.  For  the  water  dimer,  cal¬ 
culations  using  aug-pc-3,  aug-pc-2,  and  aug-pc-1  basis  sets  employ  values  of  p'  = 
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Figure  2.1:  The  water  dimer  dissociation  curve,  obtained  using  e-DFT-EE  (red, 
dot-dashed),  e-DFT-TF  (green,  dashed)  and  e-DFT-LC  (blue,  dotted).  Also  included 
are  reference  KS-DFT  results  (black,  solid),  which  are  graphically  indistinguishable 
from  the  e-DFT-EE  results.  Total  energies  are  plotted  with  respect  to  the  KS-DFT 
minimum  of  -152.430722  Hartree.  Inset,  the  curves  are  shifted  vertically  to  align  the 
energy  minima  and  horizontally  to  align  the  equilibrium  distances. 

1CT5, 10-4,  and  5  x  10-3,  respectively.  For  Li+-Be  and  CH3-CF3,  calculations  employ 
p  =  10”6.  In  each  case,  the  parameter  k  in  Eqs.  2.15  and  2.16  is  chosen  such  that 
up'  =  10  and  £  =  10-4. 


2.4.2  Water  Dimer 

Fig.  2.1  presents  the  dissociation  curve  for  the  water  dimer,  a  system  with  a  strong 
hydrogen  bond  and  significantly  overlapping  subsystem  densities.  The  curve  is  ob¬ 
tained  using  e-DFT-EE  (dot-dashed),  e-DFT-TF  (dashed),  and  e-DFT-LC  (dotted); 
KS-DFT  results  (solid)  are  also  included  for  reference.  The  e-DFT  calculations  were 
performed  using  two  embedded  subsystems,  each  corresponding  to  a  different  molecule 
in  the  dimer.  All  calculations  presented  in  the  figure  utilize  the  modified  aug-pc-3 
basis  set,  with  the  e-DFT  calculations  employing  the  supermolccular  basis  set  conven¬ 
tion.  The  dissociation  curve  is  plotted  as  a  function  of  the  oxygen-oxygen  distance, 
with  the  equilibrium  water  dimer  geometry  obtained  from  a  KS-DFT  energy  mini- 
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mization  and  with  other  geometries  obtained  by  displacing  the  two  molecules  along 
the  oxygen-oxygen  vector  while  fixing  all  other  internal  coordinates. 

The  e-DFT-EE  results  in  Fig.  2.1  agree  well  with  KS-DFT  throughout  the  range 
of  dissociation  distances.  Numerical  results  for  the  two  methods  are  graphically  in¬ 
distinguishable,  and  the  calculated  total  energies  differ  by  less  than  0.5  kcal/mol 
throughout  the  entire  attractive  branch  of  the  curve.  Exact  numerical  agreement 
between  the  e-DFT-EE  and  KS-DFT  descriptions  is  expected  only  in  the  limits  of  an 
exact  XC  functional  and  a  complete  basis  set. 

The  sensitivity  of  the  e-DFT  results  to  approximations  in  the  NAKP  is  clearly 
demonstrated  in  Fig.  2.1.  The  curve  obtained  using  e-DFT-TF  differs  significantly 
from  the  KS-DFT  reference,  exhibiting  a  dissociation  energy  that  is  underestimated 
by  40%  (~4  kcal/mol)  and  an  equilibrium  bond  length  that  is  0.15  A  too  long.  Calcu¬ 
lations  obtained  using  e-DFT-LC  are  somewhat  improved,  although  the  dissociation 
energy  is  still  overestimated  by  20%  (~2  kcal/mol)  and  the  equilibrium  bond  length 
is  underestimated  by  0.10  A.  In  the  inset  of  Fig.  2.1,  the  curvature  of  the  potential 
energy  surfaces  in  the  vicinity  of  the  minimum  are  compared,  revealing  significant  de¬ 
viations  of  the  results  obtained  using  the  approximate  NAKP  treatments  (e-DFT-TF 
and  e-DFT-LC)  with  respect  to  the  results  obtained  using  KS-DFT  and  e-DFT-EE. 

Iannuzzi  and  coworkers18  have  demonstrated  that  e-DFT  calculations  using  ap¬ 
proximate  treatments  of  the  NAKP,  including  the  TF  and  LC94  functionals,  lead  to 
qualitative  failure  in  describing  the  structure  of  liquid  water.  Fig.  2.1  illustrates  the 
origin  of  this  failure  in  terms  of  the  pairwise  interactions  among  molecules,  and  it 
suggests  that  e-DFT-EE  will  enable  the  accurate,  first-principles  simulation  of  liquid 
water  and  aqueous  solutions.  Critical  to  this  effort,  however,  is  the  efficient  and  par- 
allclizable  implementation  of  the  EE  method  for  large  systems,  which  is  discussed  in 
Section  V. 

The  sensitivity  of  the  e-DFT  calculations  to  basis  set  completeness  is  illustrated  in 
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Figure  2.2:  Basis  set  dependence  of  the  water  dimer  dissociation  curve,  illustrated 
for  calculations  using  the  (A)  modified  aug-pc-2  and  (B)  modified  aug-pc-1  basis 
sets.  Results  for  the  e-DFT-EE,  e-DFT-TF,  e-DFT-LC,  and  KS-DFT  methods  are 
reported  as  in  Fig.  2.1.  Total  energies  are  plotted  with  respect  to  the  KS-DFT 
minimum  energies  of  -152.953947  Hartree  (panel  A)  and  -152.864441  Hartree  (panel 
B). 

Fig.  2.2,  in  which  the  water  dimer  dissociation  curves  are  recalculated  using  the  mod¬ 
ified  aug-pc-2  (Fig.  2.2A)  and  modified  aug-pc-1  basis  sets  (Fig.  2.2B).  Comparison  of 
the  KS-DFT  results  and  the  e-DFT-EE  results  reveals  that  the  agreement  between  the 
methods  worsens  with  smaller  basis  set;  of  course,  both  the  KS-DFT  calculations  and 
the  e-DFT-EE  calculations  are  basis-set  dependent.  In  the  e-DFT-EE  calculations, 
smaller  basis  sets  give  rise  to  numerical  artifacts  including  the  oscillatory  behavior 
in  the  King- Handy  expression  for  the  kinetic  potential.28  For  the  modified  aug-pc-1 
basis  set  (Fig.  2.2B),  the  reasonable  agreement  between  KS-DFT  and  e-DFT-LC  is 
due  to  a  fortuitous  cancellation  of  errors  from  the  approximate  NAKP  functional  and 
small  basis  set. 


2.4.3  Li+-Be 

We  now  consider  the  heterolytic  cleavage  of  a  weak  covalent  bond,  Li+-Be— )-Li++Be, 
using  KS-DFT  and  e-DFT.  The  e-DFT  calculations  were  performed  in  the  super- 
molecular  basis  set  convention  using  two  embedded  subsystems,  one  corresponding 
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Figure  2.3:  The  Li+-Be  dissociation  curve.  Results  for  the  e-DFT-EE,  e-DFT-TF,  e- 
DFT-LC,  and  KS-DFT  methods  are  reported  as  in  Fig.  2.1.  The  results  for  e-DFT-EE 
and  the  reference  KS-DFT  results  are  graphically  indistinguishable.  Total  energies  are 
plotted  with  respect  to  the  KS-DFT  minimum  energy  of  -21.962072  Hartree.  Inset, 
the  curves  are  aligned  as  in  the  inset  of  Fig.  2.1. 

to  the  2-electron  Li  ion  and  the  other  corresponding  to  the  4-electron  Be  atom.  The 
dissociation  curve  for  Li+-Be  is  plotted  in  Fig.  2.3. 

As  is  seen  from  the  main  figure,  the  e-DFT-EE  calculations  accurately  reproduce 
the  calculated  total  energies  from  KS-DFT  throughout  the  entire  range  of  internuclear 
distances.  The  dissociation  curves  for  these  two  methods,  which  are  graphically  in¬ 
distinguishable  in  Fig.  2.3,  deviate  by  less  than  0.2  kcal/mol  throughout  the  range  of 
separations  and  the  dissociation  energy  deviates  by  only  0.07  kcal/mol.  In  contrast, 
the  e-DFT-TF  results  are  in  qualitative  disagreement  with  the  KS-DFT  reference 
calculations;  in  addition  to  dramatically  overestimating  the  dissociation  energy  of 
the  molecule  by  ~12.5  kcal/mol,  the  method  predicts  the  equilibrium  bond  length 
to  be  20%  too  short.  Interestingly,  the  e-DFT-LC  method  performs  significantly 
worse  in  this  application.  The  calculations  based  on  the  approximate  LC94  kinetic 
energy  functional  overestimate  the  dissociation  energy  by  ~16  kcal/mol  and  predict 
the  equilibrium  bond  length  to  be  25%  too  short.  The  inset  to  Fig.  2.3  illustrates 
that  both  e-DFT  methods  that  use  approximate  treatments  for  the  NAKP  lead  to 
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an  overestimation  of  the  energy  surface  curvature  in  the  vicinity  of  the  equilibrium 
bond  distance. 

The  results  in  Fig.  2.3  illustrate  the  well-known  breakdown  of  e-DFT  with  ap¬ 
proximate  treatments  of  the  NAKP  for  applications  involving  strongly  overlapping 
subsystem  densities.  They  further  show  that  our  EE  method  overcomes  this  large 
error,  yielding  the  first  numerical  demonstration  of  an  e-DFT  method  to  describe 
covalent  bond-breaking  with  chemical  accuracy.  Since  e-DFT-EE  is  a  formally  exact 
method,  this  result  is  expected.  However  demonstration  that  the  level  of  accuracy  in 
Fig.  2.3  can  be  achieved  in  practical  numerical  simulations  constitutes  a  non-trivial 
validation  of  the  method. 

2.4.4  CH3-CF3 

In  a  more  challenging  application  for  e-DFT,  we  consider  the  heterolytic  cleavage  of 
a  strong  carbon-carbon  cr-bond,  CH3-CF3  — y  CHjj"  +  CFJ.  The  e-DFT  calculations 
were  again  performed  in  the  supermolecular  basis  set  convention  using  two  embedded 
subsystems,  one  corresponding  to  the  8-electron  CH^j"  moiety  and  the  other  corre¬ 
sponding  to  the  34-electron  CF3  moiety.  The  geometry  for  the  lowest  energy  point 
along  the  curve  is  provided  in  the  supplemental  information;  the  dissociation  curve  in 
Fig.  2.4  is  plotted  by  extending  the  C-C  distance  keeping  all  other  internal  coordinates 
unchanged. 

The  dissociation  curves  in  Fig.  2.4  are  presented  only  for  e-DFT-EE  and  the  refer¬ 
ence  KS-DFT  calculations.  e-DFT-EE  reproduces  the  KS-DFT  reference  value  for  the 
total  energy  for  the  molecule  at  the  equilibrium  bond  distance  to  within  1.5  kcal/mol, 
and  the  embedding  method  also  recovers  the  reference  value  for  the  equilibrium  bond 
distance.  Furthermore,  as  is  clear  from  the  inset,  e-DFT-EE  accurately  reproduces 
the  curvature  of  the  energy  surface  in  the  vicinity  of  the  equilibrium  bond  distance.  In 
contrast,  the  e-DFT-TF  and  e-DFT-LC  descriptions  for  this  system  fail  dramatically, 
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Figure  2.4:  The  CH3-CF3  dissociation  curve  for  heterolytic  cleavage  of  the  C-C 
bond.  Results  are  presented  for  the  e-DFT-EE  (red,  dot-dashed)  and  KS-DFT  (black, 
solid)  methods.  Total  energies  are  plotted  with  respect  to  the  KS-DFT  minimum 
energy  of  -377.575687  Hartree.  Inset,  the  curves  are  aligned  as  in  the  inset  of  Fig.  2.1. 


predicting  total  energies  at  the  equilibrium  bond  distance  that  deviate  from  the  KS- 
DFT  reference  by  731  kcal/mol  and  981  kcal/mol,  respectively.  For  calculations  with 
such  strongly  interacting  subsystems,  the  failure  of  e-DFT  with  approximate  descrip¬ 
tions  for  the  NAKP  methods  has  been  previously  observed.21  However,  the  results 
for  e-DFT-EE  in  Fig.  2.4  demonstrate  significant  progress  in  the  accurate  description 
of  covalently  interacting  subsystems  using  e-DFT. 


2.5  Results:  Extension  to  Larger  Systems 

2.5.1  Pairwise  Treatment  of  the  NAKP 

In  the  previously  described  implementation  of  e-DFT-EE,  the  ZMP  step,  or  an  equiv¬ 
alent  LCS,  is  performed  on  the  full  system  of  interest.  However,  numerical  challenges 
limit  the  LCS  to  systems  with  less  than  10-15  atoms,32’33’47-50  potentially  hindering 
the  applicability  of  e-DFT-EE  in  large  systems.  To  avoid  this  problem,  we  demon¬ 
strate  an  pairwise  approximation  for  the  NAKP  that  enables  the  scalable  implemen¬ 


tation  of  e-DFT-EE. 
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For  a  system  composed  of  Nsuh  embedded  subsystems,  { pa },  the  non-additive 
kinetic  energy  can  be  approximated  using  a  pairwise  sum,24  such  that 

Nub 

U'*d[{p«}]  =  r,[p]-ET>W  Pv) 

a=l 

Nub 

~  ^  (Ts[pa  +  Pf)]  -  Ts \pa]  —  Ts[pp\) , 

a</3=l 

-Wsub 

where  p  =  The  NAKP  for  a  given  subsystem  a  is  then 

a=l 

-^sub 

^nad \Pa i  {pa}]  r]  =  Kf  (r)  -  u^(r)).  (2.18) 

Applying  the  EE  method  to  this  approximation  for  the  NAKP,  a  ZMP  step  is 
performed  at  each  iteration  of  the  KSCED  to  obtain  the  KS  orbitals  corresponding 
to  each  pair  of  subsystems  densities,  Then,  using  both  the  subsystem  KS 

orbitals  {</>"}  from  the  KSCED  and  the  subsystem-pair  KS  orbitals  the  NAKP 

is  evaluated  directly  from  Eqs.  2.8  and  2.18.  In  this  approach,  only  the  NAKP  is 
assumed  to  be  pairwise  additive;  all  other  interactions  in  the  system  are  treated 
with  full  generality.  Since  the  ZMP  step  is  applied  only  to  the  subsystem  pairs,  this 
approach  is  numerically  feasible  if  each  subsystem  is  limited  to  a  relatively  small 
number  of  atoms,  regardless  of  the  total  system  size.  The  short-ranged  nature  of 
contributions  to  the  non-additive  kinetic  energy  suggests  that  distance-based  cutoffs 
can  be  employed  within  the  sum  over  subsystem  pairs. 24 


It  was  emphasized  earlier  that  the  converged  results  of  the  ZMP  step  are  inde¬ 
pendent  of  the  choice  of  external  potential,  uext(r),  in  Eq.  2.10.  In  the  pairwise  im¬ 
plementation  of  e-DFT-EE  for  the  water  trimer  in  Sec.  V  B,  we  employ  the  following 
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external  potential  for  each  pair  of  densities  pa  and  pp, 

vext  (r)  =  vne  (r )  +  vj  [p;  r]  +  vxc  [p;  r] 

i  [p] _ STs[pa  +  pp]  ,  < 

S(pa  +  pp)  5(pa  +  pp)'  ’ 

where  Ts  indicates  the  approximate  TF  functional.  This  external  potential  approxi¬ 
mates  the  KSCED  effective  potential  (Eq.  2.3)  for  the  pair  of  subsystems  embedded 
within  the  remainder  of  the  full  system;  note  that  the  TF  functional  is  used  only  to 
regularize  the  effective  potential  for  the  ZMP  step;  it  does  not  introduce  any  addi¬ 
tional  approximation  into  the  e-DFT-EE  calculation.  In  Sec.  V  C,  we  use  a  simple 
external  potential  that  includes  only  the  electron- nuclear  interactions  for  the  subsys¬ 
tem  pair. 

The  following  two  sections  demonstrate  the  accuracy  of  this  pairwise  implemen¬ 
tation  of  e-DFT-EE  (Sec.  V  B)  and  the  efficiency  with  which  it  can  be  implemented 
in  parallel  (Sec.  V  C). 

2.5.2  Water  Trimer  Application:  Testing  Pairwise  Additivity 
in  the  NAKP 

Fig.  2.5  presents  a  test  of  pairwise  additivity  in  the  NAKP  (Eq.  2.18)  for  a  hydrogen- 
bonded  trimer  of  water  molecules.  e-DFT-EE  calculations  are  performed  using  three 
embedded  subsystems,  each  corresponding  to  a  different  molecule  in  the  trimer.  In 
a  first  set  of  results,  the  symmetric  dissociation  curve  for  the  trimer  is  calculated 
using  no  assumptions  about  the  NAKP  (solid);  in  a  second  set  of  results,  the  curve 
is  calculated  using  assuming  pairwise  additivity  of  the  NAKP  (dot-dashed).  The 
equilibrium  geometry  is  provided  in  the  supplemental  information;  other  geometries 
along  the  dissociation  curve  were  then  obtained  by  uniformly  stretching  the  oxygen- 
oxygen  distances  in  the  cluster,  keeping  all  other  internal  coordinates  unchanged. 
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Figure  2.5:  Symmetric  dissociation  curves  for  the  water  trirner,  illustrating  the 
pairwise  additivity  of  the  NAKP.  Calculations  are  performed  using  the  e-DFT-EE 
method,  with  no  approximation  to  the  NAKP  (black,  solid)  and  with  the  pairwise 
approximation  to  NAKP  (red,  dot-dashed).  The  curves  are  plotted  as  a  function  of  the 
sum  of  the  three  0-0  distances,  with  details  of  the  molecular  geometries  provided  in 
the  text.  Total  energies  plotted  with  respect  to  the  minimum  energy  of  -229.4403073 
Hartree  for  the  full  NAKP  treatment.  Inset,  the  difference  between  the  two  curves  is 
plotted. 

The  trirner  calculations  were  performed  using  the  modified  aug-pc-2  basis  set  with 
the  monomolecular  basis  set  convention;  all  other  calculation  details  are  identical  to 
those  described  previously  for  the  modified  aug-pc-2  calculations  of  the  water  dimer. 

The  agreement  between  the  two  curves  in  Fig.  2.5  indicates  that  Eqs.  2.17  and 
2.18  are  excellent  approximations  for  the  non-additive  kinetic  energy  and  NAKP, 
respectively.  Throughout  the  entire  attractive  branch  of  the  curve  the  total  energies 
differ  by  less  the  0.5  kcal/mol,  and  the  largest  deviations  appear  only  in  the  strongly 
repulsive  region  at  short  distances.  This  good  agreement  is  particularly  notable,  given 
that  the  cyclic  trirner  geometries  might  be  expected  to  magnify  possible  non-additive 
contributions  to  the  total  energy;  even  better  adherence  of  the  NAKP  to  pairwise 
additivity  is  expected  for  linear  geometries  of  the  trirner.  We  have  previously  noted 
that  higher-order  corrections  to  Eqs.  2.17  and  2.18  are  possible,24  although  the  results 
in  Fig.  2.5  suggest  that  the  assumption  of  pairwise  additivity  will  be  adequate  in  many 


cases. 
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2.5.3  Parallel  Scaling  of  e-DFT-EE 

Primary  bottlenecks  in  KS-DFT  include  calculation  of  the  two-electron  integrals  and 
solution  of  the  eigenvalue  problem.  In  standard  implementations,  the  two-electron 
integral  calculations  scales  as  M4  and  the  eigenvalue  calculation  scales  at  best  as 
M2,  where  M  is  the  total  number  of  basis  functions.51,52  More  efficient  methods  for 
computing  the  two-electron  integrals  include  prescreening, 53  Ewald  summations, 54 
and  the  fast-multipole  method; 55  however,  solution  of  the  eigenvalue  problem  remains 
a  computational  bottleneck  in  most  KS-DFT  implementations. 56 

As  has  been  noted  in  previous  work,18  the  monomolecular  basis  set  convention 
leads  to  advantageous  scaling  properties  for  e-DFT.  In  this  convention,  the  number 
of  basis  functions  used  to  solve  each  KSCED,  Msu b,  is  independent  of  system  size. 
Consequently,  the  total  cost  of  the  eigenvalue  problem  scales  linearly  with  the  number 
of  subsystems,  Nsub,  and  it  can  be  trivially  parallelized  to  the  subsystem  level. 

The  cost  of  the  two-electron  integral  calculation  is  also  reduced  in  the  monomolec¬ 
ular  basis  set  convention.  Terms  arising  from  orbitals  centered  on  molecules  in  more 
than  two  different  subsystems  are  exactly  zero,  such  that  the  total  cost  of  this  op¬ 
eration  scales  with  Ar2lbM4]b.  Furthermore,  in  this  convention,  the  density  for  each 
subsystem  is  spatially  localized,  such  that  short-ranged  contributions  to  the  KSCED 
effective  potential,  including  exchange,  correlation,  short-ranged  electrostatic  contri¬ 
butions,  and  pairwise  contributions  to  the  NAKP,  can  be  truncated  at  a  cutoff  dis¬ 
tance.  Long-ranged  electrostatic  contributions  to  the  KSCED  effective  potential  can 
be  efficiently  treated  using  Ewald  summations  or  other  standard  methods. 54,55  Setting 
aside  these  long-ranged  terms  for  the  current  demonstration,  the  use  of  distance-based 
cutoffs  reduces  the  scaling  of  the  total  two-electron  integral  calculation  to  AsubM4lb, 
which  can  be  parallelized  to  yield  constant  wall-clock  time  scaling  with  increasing 
system  size. 
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Figure  2.6:  Wall-clock  timings  for  lattices  of  hydrogen  molecules,  ranging  in  size 
from  8  to  125  H2  molecules.  The  dotted  blue  lines  indicate  ideal  quadratic  and  linear 
scaling,  the  solid,  black  curve  corresponds  to  the  serial  implementation  of  integral- 
prescreened  KS-DFT  in  Molpro,  and  the  dashed,  red  curve  corresponds  to  e-DFT-EE 
using  a  number  of  parallel  processors  equal  to  the  number  of  molecules  in  the  system. 

To  illustrate  these  scaling  properties,  Fig.  2.6  presents  benchmark  timings  for 
simple  tetragonal  lattices  of  8  to  125  H2  molecules,  using  both  e-DFT-EE  and  the 
KS-DFT  implementation  in  Molpro.  The  H2  molecules  are  oriented  parallel  to  the  z 
axis,  with  a  bond  length  of  0.8  A,  and  the  centers-of-mass  for  the  molecules  are  spaced 
by  3.0  A  along  the  x  and  y  axes  and  by  3.8  A  along  the  z  axis.  All  calculations  employ 
the  uncontracted  STO-3G  basis  set,57  Slater  exchange58  without  electron  correlation, 
and  a  grid  density  that  ensures  an  integration  error  in  the  exchange  energy  of  less 
than  10-6  Hartree.  The  e-DFT-EE  calculations  are  performed  with  each  molecule 
defined  as  a  different  subsystem,  using  the  monomolecular  basis  set  convention,  and 
using  one  parallel  processor  per  subsystem.  Values  for  the  parameters  A,  p ,  k,  and 
the  MO  shift  are  the  same  as  those  used  for  the  Li+-Be  system.  The  cutoff  for  the 
calculation  of  the  electrostatics,  exchange,  and  NAKP  terms  is  set  to  4.0  A  in  these 
calculations,  such  that  only  nearest-neighbor  molecules  in  the  lattice  contribute  to 
these  terms.  All  calculations  are  performed  on  a  cluster  of  dual,  quad-core  2.6  GHz 
Xeon  Intel  processors  with  Inhniband  communication. 

The  timings  in  Fig.  2.6  indicate  that  the  e-DFT-EE  wall-clock  time  scales  hide- 
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Figure  2.7:  Error  in  the  total  energy  of  the  e-DFT-EE  calculation  relative  to 
KS-DFT  for  increasing  system  size,  plotted  with  respect  to  the  number  of  nearest- 
neighbor  pairs. 


pendently  of  the  system  size,  with  the  deviations  at  small  sizes  due  the  boundaries  of 
the  finite  crystal.  As  expected,  the  KS-DFT  results  in  the  serial  Molpro  implementa¬ 
tion  with  integral  prescreening  scales  quadratically  with  the  increasing  system  size.  In 
Fig.  2.7,  relative  energy  of  the  e-DFT-EE  and  the  KS-DFT  calculations  are  plotted  as 
a  function  of  the  number  nearest-neighbor  pairs  in  the  lattice,  NpaiTS  =  3(A^sub  —  N^). 
The  error  is  small  and  independent  of  system  size.  The  integrated  error  in  the  density 
per  molecule  was  found  to  behave  similarly  (not  shown). 


2.6  Conclusions 

We  introduce  a  general  implementation  of  the  EE  method  for  calculating  NAKP  con¬ 
tributions  in  the  e-DFT  framework,  and  we  present  a  range  of  molecular  applications. 
The  accuracy  of  e-DFT-EE  is  demonstrated  for  systems  with  covalently  bonded  and 
hydrogen-bonded  subsystems.  For  the  dissociation  of  the  water  dimer  and  the  covalent 
bonds  in  Li+-Be  and  CH3-CF3,  e-DFT-EE  preserves  excellent  agreement  with  refer¬ 
ence  KS-DFT  calculations,  whereas  approximate  treatments  for  the  NAKP,  including 
those  based  on  the  TF  or  LC94  kinetic  energy  functionals,  lead  to  known  failures. 
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Furthermore,  pairwise  approximation  of  the  NAKP  yields  excellent  accuracy  for  the 
hydrogen-bonded  water  trimer,  and  it  enables  ideal,  constant  system-size  scaling  in 
applications  to  molecular  clusters  with  up  to  hundreds  of  atoms.  These  results  estab¬ 
lish  e-DFT-EE  as  a  promising  methodology  for  performing  accurate,  first-principles 
molecular  dynamics  and  for  accurately  embedding  high-level  wavefunction  methods 
in  complex  systems. 
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Chapter  3 

Density  functional  theory  embedding  for  correlated 
wavefunctions:  Improved  methods  for  open-shell 
systems  and  transition  metal  complexes 

3.1  Introduction 

The  demand  for  accurate  and  efficient  descriptions  of  complex  molecular  systems  re¬ 
quires  development  of  quantum  embedding  methods  for  electronic  structure  in  which 
a  small  subsystem  is  treated  with  a  high  level  of  theory  while  the  remainder  of  the 
system  is  treated  at  a  more  affordable  level.  Widely  used  examples  of  quantum 
embedding  include  QM/MM, 1-6  ONIOM,'’8  and  fragment  molecular  orbital  (FMO) 
approaches,9-11  which  have  led  to  significant  advances  in  the  simulation  of  condensed- 
phase  and  biomolecular  systems.  However,  such  methods  generally  rely  on  empirical 
models  for  the  subsystem  interactions,  including  link-atom  approximations  for  em¬ 
bedding  across  covalent  bonds12-15  and  point-charge  electrostatic  descriptions  of  the 
environment,4’9  that  are  difficult  to  systematically  improve  and  that  can  fail  in  prac¬ 
tical  applications. 3,6,16,17 

Density  functional  theory  (DFT)  offers  an  appealing  framework  for  addressing  this 
challenge. 18-42  DFT  embedding  provides  a  formulation  of  electronic  structure  theory 
in  which  subsystem  interactions  depend  only  on  their  electronic  densities,  including 
non-additive  contributions  due  to  the  electrostatic,  exchange-correlation  (XC),  and 
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kinetic  energy  terms.  In  the  WFT-in-DFT  embedding  approach,  the  DFT  embedding 
potential  is  included  as  an  external  potential  for  WFT  calculations,  providing  a  WFT- 
level  description  for  one  (or  more)  subsystem  while  the  remaining  subsystems  and 
their  interactions  are  seamlessly  treated  at  the  DFT  level  of  theory. 

Several  groups,  including  this  one,  have  recently  demonstrated  that  non-additive 

r\(?  OQ  QQ  QC 

kinetic  energy  contributions  to  the  embedding  potential  can  be  exactly  computed  ’  ’ 
with  the  use  of  optimized  effective  potential  (OEP)  methods.43  1!)  In  this  paper,  we 
introduce  a  simple  technique  to  improve  the  robustness  of  OEP  calculations  in  sys¬ 
tems  that  exhibit  small  HOMO-LUMO  gaps,  such  as  transition  metal  complexes. 

In  addition,  we  derive  spin-dependent  embedding  potentials  to  enable  the  accurate 
description  of  open-shell  systems  in  the  WFT-in-DFT  embedding  framework.  Nu¬ 
merical  applications  to  the  van-der-Waals-bound  ethylene-propylene  dimer  and  to 
the  hexaaquairon(II)  transition-metal  cation  illustrate  the  applicability  of  these  new 
techniques  and  demonstrate  the  accuracy  of  the  WFT-in-DFT  approach  in  systems 
for  which  conventional  density  functional  theory  methods  exhibit  substantial  errors. 

3.2  Theory 

Like  Kohn-Sham  (KS)-DFT,  DFT  embedding  provides  a  formally  exact  framework 
for  the  ground-state  electronic  structure  problem.  Here,  we  review  DFT-in-DFT 
embedding  and  its  basis  for  WFT-in-DFT  calculations. 

3.2.1  DFT-in-DFT  Embedding 

We  begin  by  considering  a  closed-shell  system  in  which  the  total  electronic  density 
Pab  consists  of  two  subsystems,  pab  =  Pa  +  Pb-  The  corresponding  one-electron 
orbitals  for  pa  and  pe  obey  the  Kohn-Sham  Equations  with  Constrained  Electron 


Density  (KSCED),22 
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V2  +  veS[pA,  pab;  r]  (j>f{ r)  =  eA(j)A( r),  (3.1) 

-^V2  +  veS[pB,pAB;  r]  4>f(r)  =  r),  (3.2) 

where  i  —  1, ... ,  NA ,  j  —  1, . ,  A^B,  and  NA  and  A^B  are  the  number  of  electrons  in 
the  respective  subsystems.  ves  is  the  effective  potential  for  the  coupled  one-electron 
equations, 

ves[pA,  Pab;  r]  =  p^s[pA;  r]  +  uemb(r),  (3.3) 

where  p|^s[pA;  r]  is  the  standard  KS  potential  for  subsystem  A,  and 

Vemb  (r )  =  uBc (r )  +  Vj  [pB ;  r]  +  vxc  [pAB ;  r]  - 

Dec  \pa ;  r]  +  vnad  [pA ,  pAB ;  r] .  (3.4) 

Here,  vBe(r)  is  the  nuclear-electron  Coulomb  potential  from  the  nuclei  contained  in 
subsystem  B,  vj  is  the  Hartree  potential,  and  vxc  is  the  XC  potential.  In  addition 
to  these  familiar  terms  from  conventional  KS-DFT  calculations,  DFT  embedding 
introduces  the  non-additive  kinetic  potential  (NAKP)  which  properly  enforces  Pauli 
exclusion  between  the  subsystem  densities.  It  is  obtained  from 


Diad  [pa  >  Pab  ■  P 


(3.5) 
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where  Tsnad[/?A,PB]  =  Ts[pab]  —  Ts[pa]  —  R,[pB],  and  Ts[p]  is  the  non-interacting  kinetic 
energy  functional.  The  total  energy  functional  for  the  full  system  is  then 

E[pAB]  =  Ts[pA]  +  Ts[pb]  +  rsnad[pA,  Pb]  + 

hiie  [pab]  +  J[p  ab]  +  Exc[p  ab\,  (3.6) 

where  the  last  three  terms  on  the  right-hand  side  (RHS)  are  the  nuclear-electron 
Coulomb  energy,  the  Hartree  energy,  and  the  XC  energy  computed  over  the  total 
density.  Enforcing  uemb(r)  to  be  identical  for  all  subsystems  (see  Sec.  Ill  B)  leads  to 
a  unique  partitioning  in  the  DFT  embedding  formulation,  such  that  the  specification 
of  the  nuclei  and  the  integer  number  of  electrons  in  subsystem  A  and  B  fully  specify 
the  density  partitioning. 26,35  Eqs.  1-3.6  are  easily  generalized  to  the  description  of 
multiple  embedded  subsystems. 

We  have  previously  demonstrated  that  by  using  OEP  methods  to  calculate  the 
NAKP,  DFT-in-DFT  embedding  can  accurately  describe  both  weakly  and  strongly 
interacting  subsystems,  including  subsystems  connected  by  covalent  bonds;27,28  and 
we  have  shown  that  this  method  is  computationally  feasible  for  large  systems  by  way 
of  localized  approximations  to  the  NAKP. 28  More  recently,  we  have  introduced  a 
projection  approach  that  completely  avoids  the  NAKP  calculation  in  exact  DFT  em¬ 
bedding,29  which  appears  worthy  of  further  investigation.  The  OEP-based  approach 
employed  here  is  appealing  because  it  provides  a  local  embedding  potential  that  is  a 
functional  of  only  the  subsystem  electronic  densities. 

In  practice,  the  KSCED  equations  (Eq.  3.1  and  Eq.  3.2)  are  solved  by  simply  mod¬ 
ifying  the  core  Hamiltonian  in  the  self-consistent  held  (SCF)  calculation  to  include 
the  additional  embedding  terms.  The  embedding  potential  (Eq.  3.4)  can  be  written 
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in  the  atomic  orbital  (AO)  basis  as 

Vemb  =  VB  +  J[tb]  +  Vxc[7ab]  “ 

Vxc[7a]  +  vnad[7A,  7ab],  (3.7) 

where  the  various  terms  on  the  RHS  of  this  expression  correspond  to  those  in  Eq.  3.4. 
The  subsystem  and  total  AO  density  matrices  in  Eq.  3.7  satisfy  7a  +  7b  =  7ab-  It 
follows  that  the  Fock  matrix  for  subsystem  A  can  be  expressed  as 

f A  in  B  =  hAinB  +  J[7A]+Vxc[7A])  (3.8) 


where 


hAinB  =  hA  +  vemb,  (3.9) 

and  hA  is  the  core  Hamiltonian  for  subsystem  A  (the  kinetic  energy  plus  external 
potential  due  to  the  nuclei  in  subsystem  A).  The  Fock  matrix  for  subsystem  B,  fB  m  A, 
is  analogously  defined. 

3.2.2  WFT-in-DFT  Embedding 

The  embedding  potential  in  Eq.  3.4  describes  the  subsystem  interactions  in  terms 
of  their  corresponding  electronic  densities.  However,  the  subsystem  densities  can  be 
computed  with  any  level  of  theory,  thus  allowing  for  the  description  of  one  subsystem 
at  the  (single-  or  multi-reference)  WFT  level,  while  the  remaining  environment  is 
treated  at  the  DFT  level.  23-29>34, 35, 50-54  Qosecpsiieii  WFT-in-DFT  embedding  simply 
involves  performing  a  WFT  calculation  on  a  given  subsystem  using  the  modified 
core  Hamiltonian,  hA  in  B  in  Eq.  3.9,  that  contains  the  embedding  terms  due  to  the 
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environment  of  the  other  subsystem.  The  WFT-in-DFT  energy  is  then  obtained 
by  modifying  the  DFT  energy  with  respect  to  subsystem  contributions  at  the  WFT 
level,35 


£totKFT,pEFT]  =  E»ripZT] 

-^aFT[PaFT]  +  J  vemh(r)pZT(r)dr^ 

+  (eJTIpD  +  J  ^emb(r)pAFT(r)^  •  (3-10) 

This  expression  is  easy  to  evaluate  since  the  terms  in  the  parentheses  are  just  the 
DFT  and  WFT  energies  of  subsystem  A  performed  using  the  modified  core  Hamil¬ 
tonian,  hAmB.  Just  as  DFT-in-DFT  embedding  is  an  exact  theory  for  the  case  of 
an  exact  DFT  XC  functional,  WFT-in-DFT  embedding  is  an  exact  theory  for  the 
case  of  an  exact  DFT  XC  functional  and  a  full  configuration  interaction  (FCI)  WFT 
description. 52 

3.3  Methods  of  Implementation 

Here,  we  describe  techniques  to  improve  the  accuracy  and  convergence  of  both  DFT- 
in-DFT  and  WFT-in-DFT  calculations.  First,  a  description  of  open-shell  DFT  em¬ 
bedding  is  developed  to  incorporate  the  effects  of  spin- dependence  in  the  embedding 
potential.  Then,  implementation  of  the  OEP  calculation  is  discussed,  and  an  orbital 
occupation  constraint  is  introduced  to  enable  robust  DFT-in-DFT  and  WFT-in-DFT 
embedding  calculations  for  systems  with  low-lying  virtual  orbitals,  such  as  transition 
metal  complexes. 
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3.3.1  Embedding  for  Open-Shell  Systems 

For  an  open-shell  embedded  subsystem,  the  a  and  /3  electrons  generally  experience 
different  embedding  potentials  due  to  differing  non-additive  XC  and  NAKP  con¬ 
tributions.  Previous  WFT-in-DFT  implementations  for  open-shell  systems  have  in 
practice  neglected  this  difference,  effectively  assuming  that  the  spin  polarization  is 
localized  within  the  WFT  subsystem. 34,53  In  this  study,  we  show  that  effects  due  to 
spin- dependent  potentials  are  substantial  and  easily  included  via  separate  a  and  f3 
embedding  potentials.  We  develop  approaches  to  utilize  both  restricted  and  unre¬ 
stricted  open-shell  orbital  formulations  in  WFT  calculations. 

3. 3. 1.1  Open-Shell  DFT-in-DFT  Embedding 

We  begin  by  considering  an  open-shell  system  for  which  the  total  electronic  density 
is  comprised  of  the  a  and  /3  density  of  the  two  subsystems,  /?ab  =  Pa  +  Pa  +  Pb  +  Pb- 
The  effective  potential  for  the  unrestricted  spin  orbitals2'  is 

Ves[P°Li  Pai  Pm  Pm  r]  =  Vfti’a[PAiPA,r]  +<mb(r)>  (3-11) 

where  %[S’q[Pa,  Pa!  r]  is  the  standard  KS  effective  potential  for  the  unrestricted  (U)KS 
orbitals,  and  u“mb(r)  is  a  spin-dependent  embedding  potential  applied  only  to  the  a- 
spin  electrons, 


v 


a 

emb 


^ne(r)  +  Vi  [pb;  r]  +  <[pab,  Pab;  r]  ~ 
<[Pa,  Pa;  r]  +  <ad[PA,  Pab5  r]- 


(3.12) 


The  corresponding  quantities  for  the  /3-spin  electrons  are  analogously  defined.  The 
total  energy  functional  for  the  full  open-shell  system  is  then 


s(pab]  =  r.bi]  +  r,M+iri[PA,pS]  + 

T,[pi]  +  T.IpI]  +  T^lpi  pi]  + 

Ue  [Pab]  +  ^[Pab]  +  ^xc[Pab>  PabI*  (3.13) 

Separate  OEP  calculations  are  performed  over  the  a  and  /3  spin-densities  for  the  exact 
calculation  of  the  NAKP,  which  allows  for  numerically  exact  unrestricted  open-shell 
DFT  embedding  (U-DFT-in-DFT).2' 

In  practice,  we  solve  for  the  unrestricted  spin  orbitals  by  adding  the  spin-dependent 
embedding  potentials  to  the  a  and  f3  Fock  matrices.  The  a  and  (3  embedding  potential 
can  be  written  in  the  AO  basis  as 

Vcmb  =  Vne  +  J  [7b]  +  vfc  [7aB  ,  7ab]  ~ 

VL[7a,  7a]  +  Vnad[7A,  7aBL,  (3-14) 

where  £  G  {a,  (3},  and  the  corresponding  Fock  matrices  are 

fyA  in  b  =  hA  + J[7a] +v|c[7^,74] +vfmb.  (3.15) 

The  unrestricted  spin  orbitals  for  subsystem  A  are  then  obtained  by  separately  diag¬ 
onalizing  f"'A  111  B  and  f AA  in  B  in  the  usual  way. 

Practical  implementations  for  performing  OEP  calculations  using  restricted  open- 
shell  orbitals  have  yet  to  be  developed.  We  thus  introduce  a  simple,  approximate 
scheme  for  restricted  open-shell  DFT  embedding  (RO-DFT-in-DFT).  In  this  ap¬ 
proach,  a  U-DFT-in-DFT  calculation  is  first  performed,  and  the  embedding  potentials 
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v“mb  and  vfmb  are  constructed  using  Eq.  3.14.  Then,  f“,A  in  B  and  f,:)AmB  are  con¬ 
structed  using  Eq.  3.15,  and  the  usual  RO  approach  is  employed  to  obtain  subsystem 
orbitals  that  are  spatially  identical  for  the  a  and  f3  electrons.  Specifically,  f“’A  111  B  is 
diagonalized  to  obtain  a  set  of  occupied  a  spin  orbitals,  {bAA},  and  f 1AmB  is  then 
projected  into  the  space  of  the  occupied  a  spin  orbitals  using 

f/^’A  “  B  =  cTf/3,A  in  BC;  (3.16) 

where  c  is  the  matrix  with  columns  comprised  of  the  AO  coefficients  for  {0“cc'}- 
Finally,  the  projected  Fock  matrix,  f  1,A  111  B,  is  diagonalized  to  obtain  the  set  of  RO 
orbitals,  {0ocC}>  with  the  first  N ^,A  orbitals  doubly  occupied  and  with  orbitals  A^,A  + 
1, . . . ,  Na’A  singly  occupied,  where  Na,A  and  N,8’A  indicate  the  number  of  a  and 
f3  electrons  in  subsystem  A.  Although  the  second  and  third  terms  on  the  RHS  of 
Eq.  3.15  are  updated  at  each  iteration  of  the  RO-DFT-in-DFT  calculation,  we  leave 
the  embedding  potentials  unchanged  to  avoid  performing  OEP  calculations  using 
restricted  open-shell  orbitals.  The  RO-DFT-in-DFT  energy  for  the  total  density  is 
calculated  using  Eq.  3.13. 

Several  different  schemes  have  been  proposed  to  calculate  the  embedding  poten¬ 
tial  for  open-shell  subsystems  while  neglecting  its  spin-dependence. 34,35,53  These  ap¬ 
proaches  generally  assume  that  interactions  between  the  subsystems  can  be  described 
by  a  single  embedding  potential.  For  example,  in  systems  with  an  even  number  of  elec¬ 
trons,  the  embedding  potential,  vemb  in  Eq.  3.7,  can  be  obtained  assuming  that  each 
embedded  subsystem  is  closed-shell,  and  then  the  open-shell  subsystem  is  calculated 
using  v“mb  =  Vgmb  =  Vemb-  In  this  approach,  Eq.  3.15  is  solved  self-consistently  while 
v“mb  and  Vgmb  are  held  fixed,  and  the  final  DFT-in-DFT  energy  is  calculated  using 
Eq.  3.13.  This  spin-independent  description  of  the  embedding  potential  can  be  used 
with  either  an  unrestricted  treatment  of  the  open-shell  subsystem  (U-DFT-in-DFT- 
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CS)  or  a  restricted  treatment  of  the  open-shell  subsystem  (RO-DFT-in-DFT-CS);  we 
later  employ  the  approach  to  compare  with  the  previously  described  methods  (U- 
DFT-in-DFT  and  RO-DFT-in-DFT)  that  include  spin-dependence  in  the  embedding 
potential. 

3.3. 1.2  Open-Shell  WFT-in-DFT  Embedding 

Unrestricted  open-shell  WFT-in-DFT  (U-WFT-in-DFT)  calculations  are  performed 
by  first  computing  unrestricted  Hartree-Fock  (UHF)  orbitals  in  the  spin-dependent 
embedding  potential,  and  then  using  these  orbitals  for  a  post-HF  WFT  calculation. 
The  a  and  (3  Fock  matrices  for  the  calculation  of  the  UHF  orbitals  are 

f£.A'”B  =  hA  +  J[7A]+K«[7|]+vfmb,  (3.17) 

where  £  G  {«,/?},  K  is  the  HF  exchange  matrix,  and  the  embedding  potentials 
(Eq.  3.14)  are  obtained  from  a  U-DFT-in-DFT  embedding  calculation.  The  total 
energy  is  evaluated  using 


Restricted  open-shell  WFT-in-DFT  (RO-WFT-in-DFT)  calculations  are  performed 
by  solving  for  restricted  open-shell  HF  (ROHF)  orbitals  in  the  spin-dependent  embed¬ 
ding  potential.  Just  as  in  RO-DFT-in-DFT  embedding,  a  U-DFT-in-DFT  calculation 
is  first  performed,  and  the  embedding  potentials  v“mb  and  vfmb  are  constructed  using 
Eq.  3.14.  Then,  the  Fock  matrices  in  Eq.  3.17  are  constructed  and  the  usual  approach 
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is  employed  to  obtain  RO  orbitals;  the  second  and  third  terms  on  the  RHS  of  Eq.  3.17 
are  updated  at  each  iteration  while  the  embedding  potential  is  left  unchanged.  The 
ROHF  orbitals  are  used  in  the  post-HF  WFT  calculation,  and  the  total  energy  is  then 
evaluated  using  Eq.  3.18.  We  note  that  for  the  RO-WFT-in-DFT  energy  calculation, 
the  term  EFFT[pFFT],  is  evaluated  using  the  ROKS-DFT  energy. 

3.3.2  Optimized  Effective  Potential 

As  seen  in  Eqs.  3.5  and  3.6,  DFT  embedding  requires  computation  of  both  the  kinetic 
energy,  Ts[pab],  and  its  functional  derivative.  However,  since  the  explicit  functional 
form  for  the  kinetic  energy  is  unknown,  OEP  methods  are  needed  to  obtain  these 
terms  exactly. 

The  OEP  is  the  local  potential  for  which  solution  of  the  one-electron  equations 

—  ^ V2  +  uoEp(r)  4>i  =  £i<t>i  (3.19) 

yields  orbitals  that  correspond  to  a  given  target  density  while  minimizing  the  non¬ 
interacting  kinetic  energy.  A  variety  of  methods  for  determining  such  potentials  from 
an  input  target  density  have  been  developed. 43-49  Calculations  reported  here  employ 
the  direct  optimization  procedure  developed  by  Wu  and  Yang,  3,44  in  which  the  kinetic 
energy  is  obtained  via  the  unconstrained  maximization 

Ts[pin\  =  max  {Ws  [Tdet,  uOEp(r)]}  , 

^OEp(r) 


(3.20) 
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where 

N 

2  1 

hhs  ['hdet,  VOEp(r)]  =  2  —  -  V2|0t) 

i 

+  J  (poEp(r)  -  Pm(r))ti0Ep(r)<ir 

-  CI|V«A(r)||2,  (3.21) 

and 

voEp(r)  =  ^effS[pin;  r]  +  ih(r).  (3.22) 

In  these  equations,  v\(r)  =  btgt{ r),  {<?t(r)}  comprise  an  auxiliary  basis  set  for  the 

t 

potential,  bt  are  the  corresponding  expansion  coefficients,  and  (  is  a  regularization 
parameter.  Maximization  of  Ws  utilizes  the  Newton  method  for  optimization  with 
a  back-tracking  line  search  in  the  expansion  coefficients, 55 

bb+i)  =  b(i)  +  rH-1g,  (3.23) 

where  i  is  the  iteration  number,  H  and  g  are  the  Hessian  and  gradient  of  Ws,  respec¬ 
tively,  and  r  G  [0,1]  is  the  step-size  in  the  line  search. 

In  practice,  to  obtain  the  embedding  potential,  we  do  not  explicitly  calculate  the 
NAKP  for  each  subsystem.  Instead,  for  closed  shell  subsystems,  we  directly  update 
the  embedding  potential  (Eq.  3.4)  at  each  iteration  of  the  KSCED  equations  using26 

^emb  } (r)  =  Vemb(r)  “  ^(r),  (3.24) 

where  9  G  [0,1]  is  a  damping  coefficient.  By  construction,  the  embedding  potential  for 

each  subsystem  is  identical  at  every  iteration.  The  KSCED  equations  are  initialized 
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using  i>emb(r)  =  0,  such  that  the  initial  guess  for  the  NAKP  for  subsystem  A  exactly 
cancels  the  remaining  terms  in  Eq.  3.4  (and  the  initial  guess  for  the  NAKP  for  subsys¬ 
tem  B  likewise  cancels  the  corresponding  terms).  Upon  convergence  of  the  KSCED 
equations,  uA(r)  =  0,  UggS[pAB;r]  =  u^s[pKs;r],  and  uemb(r)  =  i^b(r).  Enforcing  the 
embedding  potential  to  be  identical  for  all  subsystems  leads  to  a  unique  partitioning 
of  the  subsystem  densities. 26,35 

For  open-shell  calculations,  we  similarly  update  the  spin-dependent  embedding 
potential  (Eq.  3.12);  the  OEP  obtained  for  a  given  spin  density  is 

vOEp(r)  =  vfs’H(pl  +  Pb)>  (Pa  +  Pb);  r]  +  <4(r),  (3-25) 

and  as  in  Eq.  3.24,  uA(r)  is  used  to  update  ufmb(r)  at  each  iteration. 

Finally,  we  note  that  XC  functionals  that  include  a  fraction  of  the  exact  exchange 
can  be  employed  in  DFT  embedding  via  the  OEP  calculation.  The  HF  exchange 
matrix,  K,  is  evaluated  using  7oep,  the  OEP  density  matrix  in  the  AO  basis.  For 
DFT-in-DFT  embedding,  the  exchange  energy  is  calculated  using 

Ex  [bOEP,  7in]  =  — tr  (7oepK[7oep]) 

+tr  ((7oep  —  7m)K[7oEp]) ,  (3.26) 

where  the  second  term  on  the  RHS  corrects  the  exchange  energy  for  small  numeri¬ 
cal  differences  between  7oep  and  7;n.  For  calculations  on  the  low-spin  state  of  the 
hexaaquairon(II)  cation,  this  correction  is  found  to  be  as  large  as  20  kcal/mol;  how¬ 
ever,  the  correction  is  not  required  for  the  evaluation  of  the  WFT-in-DFT  energy 
(Eq.  3.10),  since  is  obtained  directly  from  a  KS-DFT  calculation  on  the 

full  system. 
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3.3.3  Orbital-Occupation  Freezing 

For  Ws  to  be  a  concave  function  of  uOEP(r),  it  is  necessary43  that  the  orbitals  used  to 
construct  Poep  hr  Eq.  3.21  correspond  to  the  lowest  eigenvalues  of  Eq.  3.19.  However, 
this  can  be  problematic  for  systems  with  small  energy  differences  between  the  occupied 
and  virtual  orbitals,  where  small  changes  in  VoEp(r)  can  alter  the  relative  ordering  of 
the  orbitals. 

To  illustrate  this  issue,  Fig.  3.1  shows  the  line  search  for  an  illustrative  Newton 
step  in  an  OEP  calculation  for  the  low-spin  hexaaquairon(II)  cation.  Ws  is  plotted 
as  a  function  of  r,  where  r  is  the  step-size  in  Eq.  3.23.  For  any  step-size  larger  than 
r  =  0.38  in  this  case,  the  orbital  occupancy  changes  from  one  in  which  only  t2g-like 
d  orbitals  are  occupied  to  one  in  which  eg-like  d  orbitals  are  occupied.  In  traditional 
back-tracking  line  searches,  any  step  which  increases  Ws  would  be  accepted,  including 
the  r  =  0.5  step  indicated  with  the  red  arrow.  However,  this  step  is  problematic  since 
the  Hessian  and  gradient  of  Ws  for  the  next  Newton  step  would  be  evaluated  using  a 
density  that  corresponds  to  the  wrong  orbitals.  The  net  results  are  poor  convergence 
and  incorrect  solutions  for  the  OEP. 

We  introduce  a  simple  method  to  alleviate  this  problem  by  modifying  the  back¬ 
tracking  line  search.  Reference  (r=0)  orbitals  are  computed  from  Eq.  3.19  using 
W)Ep(r)  =  vrffS[PAB;  r],  and  for  any  proposed  step-size  r,  the  corresponding  orbitals 
are  computed  using  uoEp(r)  =  [pab;  r]  +  u.\(r).  The  proposed  step  is  rejected  if 
the  overlap  between  these  two  sets  of  orbitals  is  less  than  0.5,  regardless  of  the  change 
in  Ws ;  otherwise,  it  is  subjected  to  the  usual  criteria  of  the  back-tracking  line  search 
algorithm.  Upon  rejection,  the  step-size  r  is  reduced  by  a  factor  of  2.  This  technique 
ensures  that  the  correct  orbitals  remain  occupied  throughout  the  maximization  of 
Ws .  In  Fig.  3.1,  the  proposed  step  indicated  by  the  red  arrow  is  rejected,  whereas  the 
proposed  shorter  step  indicated  by  the  black  arrow  is  accepted;  not  only  is  the  value 
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Figure  3.1:  An  illustrative  Newton  step  in  the  OEP  calculation  for  the  low- 
spin  hexaaquairon(II)  cation,  performed  with  (black)  and  without  (red)  the  orbital- 
occupation-freezing  technique.  The  technique  ensures  that  correct  orbitals  remain 
occupied  throughout  the  maximization  of  Ws.  See  text  for  details. 

of  Ws  increased,  but  the  correct  orbitals  remain  occupied.  By  utilizing  this  technique, 
we  found  that  the  maximization  of  Ws  typically  requires  less  than  20  Newton  steps 
for  the  low  spin  state  of  the  hexaaquairon(II)  cation,  whereas  the  optimization  failed 
to  converge  without  the  use  of  orbital-occupation  freezing. 

3.3.4  Computational  Details 

The  DFT  embedding  methods  employed  here  are  all  implemented  in  the  development 
version  of  the  Molpro  software  package.56  All  calculations  employ  the  supermolecular 
basis  set  convention,  in  which  the  molecular  orbitals  for  each  subsystem  are  described 
in  the  AO  basis  for  the  full  system. 57  Calculations  on  the  ethylene-propylene  dimer  use 
the  aug-cc-pVTZ  orbital  basis  set  for  the  carbon  atoms  and  the  aug-cc-pVDZ  orbital 
basis  set  for  the  hydrogen  atoms.  Calculations  on  the  hexaaquairon(II)  cation  use  the 
aug-cc-pVTZ  orbital  basis  set  for  the  iron  atom  and  the  aug-cc-pVDZ  orbital  basis 
set  for  the  hydrogen  and  oxygen  atoms.  For  the  auxiliary  basis  set  used  in  the  OEP 
calculations,  we  employ  atom-centered  Gaussian  basis  functions  (gt( r)  =  Nte~Xtr2 , 
where  Nt  is  the  normalization  constant)  for  which  the  coefficient  Xt  assumes  values 
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of  2n,  where  n  =  nmin,  nmin  +  2, . . . ,  nmax  —  2,  nmax.  Calculations  on  the  ethylene- 
propylene  dimer  employ  the  basis  set  for  which  the  s-type  functions  for  the  carbon 
and  hydrogen  atoms  span  {nmin,nmax}  =  {-4,  4},  and  the  p-type  functions  for  the 
carbon  and  hydrogen  atoms  span  {-2,  2}.  Calculations  for  the  hexaaquairon(II)  cation 
employ  the  basis  set  for  which  the  s-type  functions  for  the  iron  atom  span  {-4,  6}, 
the  p-type  functions  for  the  iron  atom  span  {-4,  6},  the  d-type  functions  for  the  iron 
atom  span  {-2,  2},  the  s-type  functions  for  the  oxygen  atoms  span  {-4,  6},  the  p-type 
functions  for  the  oxygen  atoms  span  {-2,  4},  the  s-type  functions  for  the  hydrogen 
atoms  span  {-4,  4},  and  the  p-type  functions  for  the  hydrogen  atoms  span  {-2,  2}.  For 
all  systems,  the  finite  auxiliary  basis  set  for  the  OEP  calculations  was  confirmed  to 
introduce  a  difference  of  less  than  1  kcal/mol  between  the  total  energy  computed  using 
KS-DFT  and  either  closed-shell  or  unrestricted  open-shell  DFT-in-DFT  embedding. 
The  regularization  parameter  used  in  the  OEP  calculations  is  set  to  £  =  10~3;  smaller 
values  were  tested  on  the  ethylene-propylene  dimer  and  the  hexaaquairon(II)  cation 
and  were  found  to  have  only  a  small  ((9(pHartree))  effect  on  the  total  DFT-in-DFT 
energy. 

The  KSCED  equations  are  initialized  with  subsystem  densities  comprised  of  the 
superposition  of  HF  atomic  densities  and  with  uemb(r)  =  0;  different  initial  guesses  for 
the  embedding  potential  were  tested  on  the  hexaaquairon(II)  cation  and  were  found 
to  yield  similar  final  embedding  potentials  with  only  small  ((9(10  pHartree))  changes 
in  the  total  DFT-in-DFT  energy. 
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Figure  3.2:  WFT-in-DFT  embedding  for  the  ethylene-propylene  dimer,  (a)  The 
ethylene-propylene  dissociation  curve,  obtained  using  CCSD(T)-in-B3LYP  (red)  and 
KS-DFT  with  PBE  (green),  B3LYP  (orange),  B-LYP  (blue)  and  B88-P86  (cyan)  for 
the  XC  functional.  Also  included  are  the  reference  CCSD(T)  results  (black),  which 
are  graphically  indistinguishable  from  the  CCSD(T)-in-B3LYP  results.  The  curves 
are  vertically  shifted  to  align  at  infinite  separation,  (b)  Isosurface  plots  indicate  the 
subsystem  partitioning  for  the  ethene-propene  dimer  calculations.  The  red  isosurface 
indicates  the  density  of  the  32  electrons  associated  with  the  C2H4-C2H3-  moiety,  and 
the  blue  isosurface  indicates  the  density  of  the  8  electrons  associated  with  the  -CH3 
moiety.  The  isosurface  plot  corresponds  to  an  electronic  density  of  0.05  a.u. 
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3.4  Results 

3.4.1  The  Ethylene-Propylene  Dimer:  WFT-in-DFT  Embed¬ 
ding 

The  ethylene-propylene  dimer  is  a  prototypical  system  for  which  quantum  embedding 
methods,  such  as  QM/MM  or  ONIOM,  may  be  employed.  It  exhibits  a  weak  n  —  n 
interaction  that  is  difficult  to  address  with  conventional  KS-DFT  methods,  while  also 
exhibiting  a  spectator  -CH3  moiety  that  contributes  little  to  the  interaction  energy 
while  substantially  increasing  the  cost  of  the  high-level  calculation.  However,  unlike 
the  QM/MM  treatment  of  subsystems,  the  interactions  between  the  n  —  n  system 
and  the  -CH3  moiety  can  be  treated  seamlessly  using  WFT-in-DFT  embedding,  as  is 
now  demonstrated. 

Fig.  3.2(a)  presents  the  ethylene-propylene  dimer  dissociation  curve  plotted  as 
a  function  of  the  distance  between  the  ethylene  and  propylene  tt  bonds,  with  the 
equilibrium  dimer  geometry  obtained  via  minimization  at  the  MP2  level  of  theory. 
Other  geometries  along  the  curve  are  obtained  by  displacing  the  two  molecules  along 
the  vector  formed  between  the  midpoints  of  the  two  C=C  bonds,  while  fixing  all 
other  internal  coordinates.  The  relative  energies  are  plotted  by  aligning  each  curve  at 
infinite  separation.  The  full  CCSD(T)  calculation  (black)  shows  a  binding  energy  of 
2.0  kcal/mol.  KS-DFT  calculations  using  the  PBE 58,59  (green),  B3LYP60  (orange), 
B-LYP 61,62  (blue)  and  B88-P86 61,63  (cyan)  XC  functionals  illustrate  the  difficulty  in 
describing  dispersion  interactions  using  KS-DFT.  The  PBE  functional  underestimates 
the  binding  energy  by  1.3  kcal/mol,  while  the  rest  of  the  XC  functionals  fail  to  capture 
any  of  the  attractive  interactions. 

Finally,  the  red  curve  in  Fig.  3.2(a)  presents  the  results  of  WFT-in-DFT  embed¬ 
ding,  using  a  subsystem  partitioning  in  which  the  32  electrons  associated  with  the  n 
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system  (C2H4-C2H3-,  red  in  Fig.  3.2(b))  are  treated  at  the  WFT  level  of  theory  and 
the  remaining  8  electrons  in  the  -CH3  moiety  are  treated  at  the  DFT  level  of  theory. 
We  employ  CCSD(T)  for  the  WFT  and  the  B3LYP  XC  functional  for  the  DFT  (i.e., 
CCSD(T)-in-B3LYP).  Fig.  3.2(a)  shows  excellent  agreement  between  the  CCSD(T) 
(black)  and  CCSD(T)-in-B3LYP  (red)  calculations;  these  curves,  which  are  graphi¬ 
cally  indistinguishable,  differ  by  less  than  0.10  kcal/mol  through  the  entire  range  of 
distances.  We  have  confirmed  that  this  level  of  accuracy  is  maintained  with  different 
XC  functionals  used  for  the  DFT;  specifically,  CCSD(T)-in-(B-LYP)  energies  differ 
from  the  CCSD(T)  results  by  less  than  0.20  kcal/mol  throughout  the  entire  curve. 
These  results  illustrate  that  WFT-in-DFT  embedding  can  be  used  to  systematically 
improve  DFT  results  and  to  avoid  embedding  errors  while  partitioning  across  covalent 
bonds. 

3.4.2  The  Hexaaquairon(II)  Cation 

We  now  present  DFT-in-DFT  and  WFT-in-DFT  calculations  for  the  high-spin  [5T2g  : 
(t2g)4(eg)2]  and  low-spin  [1Alg  :  (t2g)6(eg)0]  states  of  the  hexaaquairon(II)  cation,  a 
system  that  presents  challenges  due  to  the  presence  of  low-lying  unoccupied  orbitals, 
the  important  role  of  unpaired  electrons,  and  the  relatively  large  number  of  electrons 
(84  e“)  in  the  full  system.  First,  we  test  the  accuracy  of  DFT-in-DFT  embedding  for 
the  various  treatments  of  the  open-shell  embedding  potential  described  earlier.  We 
then  employ  WFT-in-DFT  calculations  to  investigate  the  low-spin/high-spin  energy 
splitting  and  the  ligation  energy  for  this  transition  metal  complex. 

3.4. 2.1  DFT-in-DFT  Embedding 

Fig.  3.3(a)  presents  the  potential  energy  curve  for  the  simultaneous  dissociation  of  all 
six  H2O  ligands  of  the  hexaaquairon(ll)  cation,  plotted  as  a  function  of  the  average 
iron-oxygen  distance.  The  equilibrium  geometries  for  the  low-spin  [1Alg  :  (t2g)6(eg)0] 
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Figure  3.3:  DFT-in-DFT  embedding  for  the  hexaaquairon(II)  cation,  (a)  The  po¬ 
tential  energy  curve  for  the  simultaneous  dissociation  of  the  six  H2O  ligands.  All 
curves  in  the  main  panel  are  vertically  shifted  to  share  a  common  minimum  energy; 
they  are  not  horizontally  shifted.  The  dissociation  curves  for  the  low-spin  (1Aig)  state 
obtained  using  KS-DFT  (black)  and  DFT-in-DFT  (red)  are  graphically  indistinguish¬ 
able.  The  dissociation  curves  or  the  high-spin  (5T2g)  state  obtained  using  UKS-DFT 
(blue),  U-DFT-in-DFT  (green),  ROKS-DFT  (magenta),  and  RO-DFT-in-DFT  (or¬ 
ange)  are  likewise  graphically  indistinguishable.  The  inset  shows  these  four  high-spin 
potential  energy  curves,  with  each  curve  vertically  shifted  only  by  the  UKS-DFT  min¬ 
imum  energy  of  —1721.693423  Hartree.  The  dashed  black  dissociation  curve  in  the 
main  panel  is  obtained  using  the  RO-DFT-in-DFT-CS  method,  which  neglects  spin- 
dependence  in  the  embedding  potential,  (b)  Isosurface  plots  indicate  the  subsystem 
partitioning  for  the  hexaaquairon(II)  cation.  The  red  isosurface  indicates  the  density 
of  the  24  electrons  associated  with  the  Fe  atom,  and  the  blue  isosurface  indicates  the 
density  of  the  60  electrons  associated  with  the  six  H20  ligands.  The  isosurface  plot 
corresponds  to  an  electronic  density  of  0.05  a.u. 


and  high-spin  [5T2g  :  (t2g)4(eg)2]  states  are  obtained  using  KS-DFT  energy  minimiza¬ 
tion  with  the  B3LYP  XC  functional;  all  other  geometries  are  obtained  by  uniformly 
stretching  the  iron-oxygen  distances  in  the  complex,  keeping  all  other  internal  co- 
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ordinates  unchanged.  All  KS-DFT  and  DFT-in-DFT  embedding  results  reported  in 
this  section  are  obtained  using  the  B3LYP  XC  functional.  The  curves  in  the  main 
panel  of  Fig.  3.3(a)  are  vertically  shifted  to  share  a  common  minimum  value;  they  are 
not  horizontally  shifted.  The  high-spin  state  is  lower  in  energy  and  exhibits  a  longer 
average  iron-oxygen  distance  than  the  low-spin  state. 

We  perform  DFT-in-DFT  embedding  using  a  subsystem  partitioning  in  which 
the  24  electrons  associated  with  the  iron  center  comprise  one  subsystem  (red  in 
Fig.  3.3(b))  and  the  remaining  60  electrons  associated  with  the  six  water  ligands 
comprise  a  second  subsystem  (blue  in  Fig.  3.3(b)).  For  the  low-spin  state,  Fig.  3.3(a) 
demonstrates  good  numerical  agreement  between  DFT-in-DFT  (red)  and  KS-DFT 
(black);  the  relative  energies  differ  by  less  than  0.6  kcal/mol  throughout  the  range  of 
reported  internuclear  distances. 

For  the  high-spin  state  of  the  hexaaquairon(II)  cation,  Fig.  3.3(a)  shows  that  the 
UKS-DFT  and  ROKS-DFT  methods  are  in  good  agreement  with  each  other,  as  well  as 
with  the  corresponding  U-DFT-in-DFT  and  RO-DFT-in-DFT  embedding  approaches 
described  in  Sec.  Ill  A  1.  The  U-DFT-in-DFT  calculation  accurately  reproduces  the 
relative  energies  obtained  from  UKS-DFT  to  within  0.4  kcal/mol  throughout  the 
attractive  branch  of  the  curve  and  to  within  0.8  kcal/mol  at  shorter  distances.  The 
RO-DFT-in-DFT  calculation  reproduces  the  relative  energy  obtained  from  ROKS- 
DFT  to  within  1.0  kcal/mol  throughout  the  attractive  branch  of  the  curve  and  to 
within  2.2  kcal/mol  at  shorter  distances. 

The  inset  of  Fig.  3.3(a)  shows  the  various  potential  energy  curves  computed  for 
the  high-spin  state  of  the  hexaaquairon(II)  cation,  with  each  curve  vertically  shifted 
by  only  the  UKS-DFT  minimum  energy.  This  inset  demonstrates  relatively  small 
differences  in  the  total  energies  computed  with  the  various  embedding  and  open-shell 
treatments. 

Finally,  the  dashed  black  curve  in  Fig.  3.3(a)  demonstrates  the  importance  of 
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including  spin-dependence  in  the  embedding  potential.  This  curve  corresponds  to  the 
RO-DFT-in-DFT-CS  treatment  of  the  high-spin  state  of  the  hexaaquairon(II)  cation 
described  in  Sec.  Ill  A  1.  It  exhibits  large  relative  errors  (over  70  kcal/mol)  compared 
to  the  other  treatments  of  the  high-spin  state  of  the  hexaaquairon(II)  cation,  as 
well  as  qualitatively  incorrectly  shortening  of  the  equilibrium  internuclear  distance. 
Although  this  approximation  is  expected  to  be  more  reliable  for  systems  in  which  the 
spin-density  is  strongly  localized  with  a  single  subsystem,  the  result  demonstrates  that 
substantial  errors  can  emerge  clue  to  the  neglect  of  spin-dependence  in  the  embedding 
potential. 

3. 4. 2. 2  WFT-in-DFT  Embedding 

We  now  consider  WFT-in-DFT  embedding  for  the  hexaaquairon(II)  cation,  employ¬ 
ing  the  same  subsystem  partitioning  as  in  the  DFT-in-DFT  embedding  calculations 
(Fig.  3.3(b)).  The  hexaaquairon(II)  cation  is  a  benchmark  system  for  spin  splittings 
in  transition  metal  complexes.66  We  initially  discuss  results  for  MP2  embedding  to 
compare  the  U-WFT-in-DFT  and  RO-WFT-in-DFT  approaches,  and  we  then  present 
results  obtained  using  CCSD(T)  embedding. 

Fig.  4.2  presents  results  for  the  low-spin/high-spin  energy  difference  (AElh)  ob¬ 
tained  using  MP2,  KS-DFT,  and  MP2-in-DFT  embedding;  detailed  values  are  re¬ 
ported  in  Table  3.1.  For  KS-DFT  calculations  of  AElh,  the  energy  for  the  high-spin 
state  of  the  hexaaquairon(II)  cation  was  obtained  at  the  UKS-DFT  level  of  theory. 
The  WFT-in-DFT  embedding  energy  for  the  low-spin  state  of  the  hexaaquairon(II) 
cation  is  obtained  using  closed-shell  WFT-in-DFT  (Sec.  II  B),  while  the  high-spin 
state  is  treated  using  either  U-WFT-in-DFT  or  RO-WFT-in-DFT  (Sec.  Ill  A  2).  The 
KS-DFT  results  (red  in  Fig.  4.2)  exhibit  strong  dependence  on  the  XC  functional, 
with  hybrid  functionals  underestimating  AElh  to  a  somewhat  lesser  degree  than  the 


semi-local  functionals. 
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Fig.  4.2  clearly  illustrates  that  the  R0-MP2-in-DFT  results  (blue)  are  in  better 
agreement  with  the  full  MP2  calculation  than  the  corresponding  U-MP2-in-DFT  re¬ 
sults  (green),  particularly  for  semi-local  XC  functionals.  Removal  of  spin-contamination 
in  the  WFT  calculation  reduces  the  energy  of  the  high-spin  state  RO-WFT-in-DFT 
calculation  with  respect  to  that  obtained  using  U-WFT-in-DFT. 

Another  important  observation  from  Fig.  4.2  is  that  the  dependence  of  AElh  on 
the  DFT  XC  functional  is  greatly  reduced  in  the  embedding  calculation,  even  though 
only  the  single  transition  metal  atom  is  treated  at  the  WFT  level.  The  spread  of 
values  obtained  at  the  KS-DFT  level  of  theory  is  over  6000  cm-1,  which  is  reduced 
by  a  factor  of  3  in  the  RO-MP2-in-DFT  embedding  calculations. 


Figure  3.4:  MP2-in-DFT  embedding  for  the  hexaaquairon(II)  cation.  High- 
spin/low-spin  splitting  energies  obtained  using  KS-DFT  (red,  circles),  U-MP2-in-DFT 
(green,  squares),  and  RO-MP2-in-DFT  (blue,  triangles)  with  a  range  of  different  XC 
functionals  that  include  B-LYP, 61,62  PBE, 58,59  PW91,64  B3LYP,60  and  PBEO.65  The 
black  line  indicates  the  reference  value  of  16439  cm-1  obtained  at  the  RO-MP2  level 
of  theory;  U-MP2  yields  a  value  of  17396  cm-1. 


Fig.  3.5(a)  presents  calculations  of  the  low-spin/high-spin  splitting  obtained  us¬ 
ing  WFT-in-DFT  calculations  at  the  RO-CCSD(T)-in-DFT  level  of  theory;  detailed 
values  are  reported  in  Table  3.2.  For  the  reference  calculation  obtained  at  the  full 
RO-CCSD(T)  level  of  theory,6'  no  T2  amplitudes  were  found  to  exceed  0.05,  indicat¬ 
ing  that  a  single-reference  description  of  the  wavefunction  is  adequate.  The  general 
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Table  3.1:  High-spin/low-spin  splitting  energies  in  cm-1  for  the  hexaaquairon(II) 
cation  obtained  using  KS-DFT,  U-MP2-in-DFT,  and  RO-MP2-in-DFT  with  a  range 
of  different  XC  functionals. 


Functional 

KS-DFT 

U-MP2-in-DFT 

RO-MP2-in-DFT 

B-LYP 

7828 

12604 

15294 

PBE 

9479 

11079 

13395 

PW91 

8593 

10924 

13201 

B3LYP 

11206 

14387 

14703 

PBEO 

14154 

13812 

13979 

RO-MP2  yields  16439  cm  1 
U-MP2  yields  17396  cm'1 


trend  for  the  RO-CCSD(T)-in-DFT  calculations  is  consistent  with  the  results  ob¬ 
tained  from  R0-MP2-in-DFT.  It  is  again  seen  that  the  dependence  of  ALfm  on  the 
XC  functional  is  substantially  reduced  using  RO-CCSD(T)-in-DFT  embedding,  and 
the  accuracy  of  the  KS-DFT  results  are  generally  improved  by  treating  the  transition 
metal  atom  at  the  WFT  level.  For  this  system,  the  embedded  RO-CCSD(T)  cal¬ 
culation  involves  correlating  significantly  fewer  electrons  than  the  full  RO-CCSD(T) 
calculation,  and  we  found  that  the  WFT  step  in  the  RO-CCSD(T)-in-DFT  calcula¬ 
tion  required  approximately  50  times  less  wall-clock  time  than  the  full  RO-CCSD(T) 
calculation. 

Fig.  3.5(b)  shows  that  the  LDA  functional68,69  presents  an  interesting  outlier  com¬ 
pared  to  the  other  results  in  Fig.  3.5(a).  Unlike  the  semi-local  and  hybrid  functionals, 
RO-CCSD(T)-in-LDA  calculations  do  not  exhibit  a  significant  improvement  with  re¬ 
spect  to  the  corresponding  KS-DFT  result.  We  now  show  that  this  anomalous  result 
arises  from  a  density-based  error  in  the  LDA  functional. 

Fig.  3.5(c)  and  Fig.  3.5(d)  present  the  charge  on  the  Fe  atom  from  a  Mulliken 
population  analysis  for  the  low-spin  state  of  the  hexaaquairon(II)  cation.  Fig.  3.5(c) 
shows  that  the  semi-local  and  hybrid  functionals  all  yield  a  similar  charge  for  the  Fe 
atom,  which  is  very  close  to  that  of  the  full  (relaxed)  CCSD  density.  In  contrast, 
Fig.  3.5(d)  reveals  the  LDA  functional  significantly  underestimates  the  Fe  atomic 
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charge,  which  indicates  a  significant  error  in  the  calculation  of  the  ground  state  den¬ 
sity.  Although  the  use  of  embedded  WFT  can  be  expected  to  overcome  the  error  in 
the  contribution  to  the  spin-splitting  energy  due  to  the  LDA  functional,  it  cannot 
overcome  this  error  in  the  actual  ground  state  density  due  to  LDA. 


B-LYP  PBE  PW91  B3LYP  PBEO  LDA  B-LYP 


Figure  3.5:  CCSD(T)-in-DFT  embedding  for  the  hexaaquairon(II)  cation.  (a,b) 
High-spin/low-spin  splitting  energies  obtained  using  KS-DFT  (red,  circles)  and  RO- 
CCSD(T)-in-DFT  (blue,  triangles)  with  a  range  of  different  XC  functionals.  The 
B-LYP+LDA  result  is  obtained  using  the  B-LYP  XC  functional  for  the  density  cal¬ 
culation  and  the  LDA  XC  functional  for  the  energy  calculation,  as  is  described  in 
the  text.  The  black  line  indicates  the  reference  value  of  14149  cm-1  obtained  at  the 
RO-CCSD(T)  level  of  theory.  (c,d)  The  charge  on  the  Fe  atom  is  obtained  using  the 
Mulliken  population  analysis  of  the  KS-DFT  calculation  with  each  functional.  The 
relaxed  CCSD  density,  indicated  by  the  black  line,  has  an  Fe  atomic  charge  of  2.56. 


To  confirm  this  interpretation,  we  show  that  removing  the  error  in  the  LDA  den¬ 
sity  leads  to  improved  WFT-in-DFT  estimates  for  the  spin-splitting  energy,  even  if 
the  LDA  functional  is  still  employed  for  the  DFT  contributions  to  the  energy.  In 
Fig.  3.5(b),  the  B-LYP+LDA  result  for  WFT-in-DFT  embedding  (blue,  triangle)  is 
obtained  by  (i)  calculating  the  embedding  potential  and  the  subsystem  densities  us- 


Table  3.2:  High-spin/low-spin  splitting  energies  in  cm-1  for  the  hexaaquairon(II) 
cation  obtained  using  KS-DFT  and  RO-CCSD(T)-in-DFT  with  a  range  of  different 
XC  functionals. 


Functional 

KS-DFT 

RO-CCSD(T)-in-DFT 

B-LYP 

7828 

12554 

PBE 

9479 

11238 

PW91 

8593 

10712 

B3LYP 

11206 

12634 

PBEO 

14154 

12912 

RO-CCSD(T)  yields  14149  cm'1. 


ing  the  B-LYP  XC  functional,  (ii)  performing  the  embedded  WFT  calculation  at  the 
CCSD(T)  level,  and  (in)  using  the  LDA  functional  and  CCSD(T)  to  evaluate  the 
respective  DFT  and  WFT  contributions  to  the  total  energy  in  Eq.  3.18.  The  cor¬ 
responding  B-LYP+LDA  result  for  KS-DFT  (red,  circle)  is  obtained  by  calculating 
the  total  density  using  KS-DFT  with  the  B-LYP  XC  functional  and  then  using  the 
LDA  functional  to  evaluate  the  KS-DFT  energy.  As  is  seen  in  Fig.  3.5(d),  the  B-LYP 
treatment  of  the  subsystem  densities  leads  to  the  expected  partial  charge  for  the  Fe 
atom;  it  avoids  the  error  in  the  electronic  density  that  is  introduced  using  LDA.  How¬ 
ever,  the  spin-splitting  energy  obtained  using  the  B-LYP+LDA  result  for  KS-DFT 
is  essentially  no  better  than  that  obtained  using  KS-DFT  with  the  LDA  functional 
(Fig.3.5(b)),  indicating  that  simply  correcting  the  LDA  error  in  the  density  is  not 
enough  to  avoid  the  LDA  error  in  the  energies.  Finally,  Fig.  3.5(b)  shows  that  the 
B-LYP+LDA  result  for  WFT-in-DFT  does  exhibit  a  substantial  improvement  over 
the  corresponding  KS-DFT  result;  this  confirms  that  WFT  embedding  is  able  to  over¬ 
come  energy-based  errors  due  to  the  DFT  XC  functional,  although  it  is  less  effective 
at  overcoming  density-based  errors  due  to  the  DFT  XC  functional. 

Although  we  have  shown  that  WFT-in-DFT  embedding  with  the  subsystem  par¬ 
titioning  shown  in  Fig.  3.3(b)  generally  leads  to  improved  estimates  for  the  low- 
spin/high-spin  splitting  energy  over  KS-DFT,  the  same  does  not  hold  true  for  calcu- 
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lated  ligation  energies  of  the  hexaaquairon(II)  cation.  Ligation  energies  calculated  us¬ 
ing  RO-CCSD(T)-in-DFT  embedding  are  essentially  unchanged  from  those  obtained 
using  KS-DFT  with  the  corresponding  XC  functional;  indeed,  the  mean  absolute  dif¬ 
ference  between  the  computed  WFT-in-DFT  and  KS-DFT  ligation  energy  is  only  0.6 
kcal/mol  per  ligand  across  the  set  of  functionals  that  includes  LDA,  B-LYP,  PBE, 
PW91,  B3LYP,  and  PBEO.  Unlike  the  spin-splitting  energy,  which  is  highly  sensitive 
to  the  electronic  structure  of  the  Fe  atom  and  is  thus  impacted  by  the  WFT  sub¬ 
system  description,  the  ligation  energy  is  dominated  by  interactions  between  the  Fe 
atom  and  the  water  ligands;  these  inter-subsystem  interactions  are  still  treated  es¬ 
sentially  at  the  DFT  level  in  WFT-in-DFT  embedding.  An  improved  description  for 
the  ligation  energy  could  be  obtained  by  simply  expanding  the  number  of  electrons 
that  are  treated  at  the  WFT  level  of  theory,  or  by  including  two-body  correlation 
corrections  through  an  embedded  many-body  expansion  description  of  the  system. 29 

3.5  Conclusion 

In  this  work,  we  have  introduced  and  demonstrated  improved  methods  for  the  imple¬ 
mentation  of  WFT-in-DFT  calculations  for  open-shell  systems  and  systems  with  low- 
lying  virtual  orbitals.  A  simple  orbital-occupation-freezing  technique  is  introduced  to 
enable  robust  OEP  calculations  on  systems  with  small  HOMO-LLIMO  gaps,  leading  to 
accurate  DFT-in-DFT  and  WFT-in-DFT  embedding  calculations  on  transition-metal 
complexes.  Furthermore,  the  use  of  spin-dependent  embedding  potentials  is  shown  to 
preserve  the  accuracy  of  open-shell  DFT-in-DFT  calculations  in  both  the  restricted 
and  unrestricted  orbital  formulations,  whereas  neglect  of  the  spin  polarization  leads 
to  significant  errors  in  both  computed  energies  and  geometries.  WFT-in-DFT  calcu¬ 
lations  on  the  hexaaquairon(II)  cation  reveal  that  the  treatment  of  only  the  single 
transition  metal  atom  leads  to  significant  improvements  in  the  accuracy  of  calculated 


spin-splittings,  as  well  as  marked  reduction  in  the  dependence  of  results  on  the  DFT 
XC  functional.  Taken  together,  the  exact  embedding  techniques  reported  and  demon¬ 
strated  here  offer  a  promising  approach  to  the  robust  treatment  of  systems  for  which 
the  accuracy  of  WFT  is  required  but  for  which  the  cost  of  the  full  WFT  calculation 
is  not  feasible. 
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Chapter  4 

Accurate  and  systematically  improvable  density  func¬ 
tional  theory  embedding  for  correlated  wavefunc- 
tions 

4.1  Introduction 

The  observation  that  many  chemical  processes  are  predominately  governed  by  changes 
within  a  localized  subsystem  has  motivated  the  development  of  a  number  of  multi¬ 
scale  strategies. 1-21  The  success  of  such  methods  is  contingent  on  the  availability  of 
a  sufficiently  accurate  description  of  the  environment,  as  well  as  a  suitable  model 
for  the  coupling  between  subsystems.  Density  functional  theory  (DFT)  provides  an 
ideal  framework  for  multiscale  embedding. 1  24  In  these  approaches,  an  electronic 
structure  calculation  on  a  chemical  system  is  partitioned  into  calculations  on  two 
subsystems:  subsystem  A,  which  is  treated  using  an  accurate  wavefunction  theory 
(WFT),  and  subsystem  B,  which  is  treated  using  the  more  computationally  efficient 
DFT  method. 25-36  Our  projector-based  WFT-in-DFT  embedding  approach  has  the 
advantage  of  offering  a  framework  that  is  both  exact  for  cases  in  which  both  subsys¬ 
tems  are  treated  using  DFT  (DFT-in-DFT  embedding)  and  efficient  for  calculations 
on  large  systems. 36,37 

Although  projector-based  embedding  is  numerically  exact  for  DFT-in-DFT  em¬ 
bedding,  it  is  clear  that  some  error  is  introduced  into  any  practical  WFT-in-DFT 
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embedding  calculation.  Because  the  energy  of  the  DFT  environment  is  calculated  at 
the  DFT  level,  this  contribution  will  be  no  more  accurate  than  that  of  a  standard  DFT 
calculation.  Evaluation  of  the  interaction  between  subsystems  is  also  handled  using 
DFT  theory,  which  introduces  errors  into  both  the  embedding  potential  of  the  WFT 
subsystem,  and  the  nonadditive  energy  between  subsystems.  We  analyze  WFT-in- 
DFT  embedding  by  decomposing  the  error  into  these  three  contributions,  and  use  the 
results  to  suggest  further  improvements  to  projector-based  embedding.  The  analysis 
is  performed  through  careful  comparison  with  local  coupled-cluster  calculations. 

We  also  analyze  the  errors  of  a  number  of  embedding  calculations  on  systems  that 
might  be  expected  to  be  particularly  difficult  to  treat  using  projector-based  embed¬ 
ding.  In  particular,  we  investigate  the  potential  energy  surface  of  a  heterolytic  bond 
cleavage  using  projector-based  embedding.  As  with  other  local  correlation  methods, 
our  embedding  method  exhibits  discontinuities  in  the  potential  energy  surface;  how¬ 
ever,  these  discontinuities  are  small  and  decrease  as  the  WFT  subsystem  is  expanded. 
Finally,  we  consider  reactions  involving  highly  conjugated  molecules,  and  find  that 
projector-based  embedding  produces  reliably  accurate  results  for  reactions  involving 
moderate  changes  in  polarization. 

4.2  Projector-Based  Embedding 

The  projector-based  embedding  method  provides  a  rigorous  framework  for  embedding 
a  DFT  or  WFT  subsystem  description  in  a  DFT  environment.36  In  this  approach, 
a  Kohn-Sham  (KS)-DFT  calculation  is  first  performed  over  the  full  system.  The 
resulting  occupied  molecular  orbitals  (MOs),  {</>*},  are  then  localized  and  partitioned 
into  the  sets  {0A}  and  {0 f },  which  correspond  to  subsystems  A  and  B,  respectively. 
These  two  sets  of  orbitals  are  used  to  form  the  density  matrices  of  subsystems  A  and 
B  in  the  atomic  orbital  basis,  yA  and  yB. 
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Next,  the  subsystem  Fock  matrix  is  formed  for  the  embedding  calculation,  such 
that 


fA  =  hAinB[7A;7B]+g[~AL  (4.1) 

where  the  embedded  core  Hamiltonian  is 


hA  111  B[yA,  7B]  =  h  +  g[7A  +  7B]  -  g[7A]  +  pPB.  (4.2) 

Here,  h  is  the  standard  one-electron  core  Hamiltonian,  g  includes  all  the  two-electron 
terms.  PB  is  a  projection  operator,  and  ji  is  a  level-shift  parameter;  7A  is  the  density 
matrix  associated  with  the  MO  eigenstates  of  fA,  {0A}.  The  projection  operator  is 
given  by 

El^fl  IW,  (4.3) 

l  ieB  J 

where  a,  (3  label  the  atomic  orbital  basis  functions.38-  5  In  the  limit  of  /i  — »  oo,  the 
MOs  in  {0A}  are  constrained  to  be  mutually  orthogonal  with  the  MOs  of  subsystem 
B;36,3'  if  in  addition  the  same  density  functional  is  used  for  all  calculations,  the  MOs 
{0A}  coincide  with  the  original  orbitals  {</>A}. 

A  self-consistent  field  optimization,  using  the  Fock  matrix  fA,  is  performed  to 
obtain  yA,  and  the  final  DFT-in-DFT  energy  is 

£dft[7A;7A,7B]  = 

^dft[7A]  +  -Soft  [7B]  +  -^dft[7A>  7B]  (4-4) 

+  tr  [(7A-7A)(hAinB[7A,7B]-h)], 

where  Pdft  is  the  standard  DFT  energy  (evaluated  with  core-Hamiltonian  h)  and 
-^dft[7A>  7B]  is  fiie  nonadditive  energy  between  the  subsystem  densities.  The  last 
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term  is  a  first-order  correction  to  the  difference  between  -Ej)^[7A,  7s]  and  77£,ft[7A>  7B] 
In  the  limit  of  /j,  — >  00,  7A  =  7A  and  the  DFT-in-DFT  embedding  energy  is  identical 
to  the  energy  from  the  corresponding  KS  calculation  performed  over  the  full  system; 
as  a  result,  the  projector-based  approach  is  numerically  exact  for  DFT-in-DFT  em¬ 
bedding  calculations.36  In  practice,  a  large  finite  value  of  p  is  used,  and  an  additional 
perturbative  correction  to  the  energy  can  performed;36  for  appropriate  values  of  p 
this  correction  is  typically  far  smaller  than  the  energy  differences  discussed  in  this 
paper  and  is  thus  neglected  throughout.  Furthermore,  as  has  been  previously  em¬ 
phasized,  this  embedding  scheme  is  exact  for  any  self-consistent  held  method,  such 
as  Hartree-Fock  (HF)  theory. 36,37 

The  nonadditive  contribution  to  the  energy,  -Ej^T[7A,  yB],  can  be  decomposed 
into  electrostatic  and  exchange-correlation  contributions 

S£ft[tA.  7B]  =  -/"“V,  7B]  +  JCV,  7B],  <4-5) 

where 


and 

£xcd[7A,  7B]  =  Ex c[7a  +  7B]  -  £xc[7a]  -  £xc[7b].  (4.7) 

The  electrostatic  term,  Jnad,  is  easily  evaluated,  and  although  the  exact  form  of  Exc 
is  unknown,  approximate  functionals  are  well  established.  Since  the  embedded  MOs 
{0A}  are  orthogonal  to  those  in  subsystem  B,  there  is  no  nonadditivity  in  the  kinetic 
energy.  This  removes  the  requirement  of  performing  optimized  effective  potential 
calculations 20,21,23,2  ,30,31  or  using  approximate  nonadditive  kinetic  energy  functionals. 
The  projector-based  formalism  easily  allows  for  WFT-in-DFT  embedding,  in  which 


Jn ad[7A,7B]=  /  dr,  /  dr; 


7a(1)7b(2) 
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(4.6) 


31 


subsystem  A  is  treated  using  a  WFT-level  description  and  subsystem  B  is  described  at 
the  DFT  level.36  This  is  achieved  by  replacing  the  standard  one-electron  core  Hamil¬ 
tonian  with  the  embedded  core  Hamiltonian  of  Eq.  4.2.  The  electronic  energy  from 
the  WFT-in-DFT  embedding  approach  is 

^wft[^a;7a,7b]  =  (TA|ilAinB[7A,7B]|TA> 

-tr[7A(hAinB[7A,7B]-h)] 

(48) 

+  -Soft  [7B] 

+  £-dT[7A,7B], 

where  |TA)  is  the  embedded  wavefunction  from  the  WFT  method,  and  HA  111  B[7A,  7B] 
is  the  WFT  Hamiltonian  resulting  from  replacing  the  standard  one-electron  core 
Hamiltonian  with  the  embedded  core  Hamiltonian.  The  term  tr  [7A(hA  111  B[7A,  7B]  —  h 
is  included  in  the  first  term  of  Eq.  4.8  and  thus  does  not  show  up  in  the  first-order 
correction  term,  as  it  did  in  Eq.  4.4. 

4.3  Results  I:  Sources  of  error  in  WFT-in-DFT 
embedding 

4.3.1  Term-By-Term  Comparison  with  LCSSD(T) 

We  now  formulate  an  approach  to  compare  the  individual  terms  in  the  energy  ex¬ 
pression  of  a  CCSD(T)-in-DFT  embedding  calculation  with  the  corresponding  val¬ 
ues  calculated  at  the  CCSD(T)  level.46  To  do  this,  we  first  recognize  that  the  local 
(L)CCSD(T)  method  by  Schutz  and  Werner 47,50-52  becomes  exactly  equivalent  to  the 
canonical  CCSD(T)  method  when  all  orbital  pairs  are  correlated  and  all  excitation 
domains  are  set  to  the  full  virtual  basis.  The  terms  in  the  LCCSD(T)  energy  expres- 
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sion,  in  turn,  can  be  organized  in  a  way  that  enables  direct  comparison  to  the  terms 
in  the  CCSD(T)-in-DFT  embedding  energy  expression. 

The  LCCSD(T)  energy  can  be  decomposed  as  a  function  of  the  amplitudes  and 
the  atomic-orbital  density  matrices  as 


^lccsd(t)[Ti,T2,T3;7A,7B]  — 

+  -E'HF  [bA]  +  £fs)  +  ^(D)  +  T) 

(4.9) 

+  £Hf[7B]  +  ^fs)  +  ^(D)  +  £(t) 

+  ^Fd[7A,  7B]  +  £(D?  +  E(T)1 

where  Eu p  is  the  HF  energy  and  EmFd[yA,7B]  is  the  same  as  Eq.  4.5,  except  with  the 
corresponding  exchange  terms  replacing  E°)ld[7A,  yB].  When  the  full  virtual  space  is 
included,  the  singles  are  additive  and  thus  there  is  no  nonadditive  component.  The 
nonadditive  correlation  for  the  double-excitation  terms  is  simply 

qS)  =  Em  -  77  -  b(bD)  (4.io) 

and  likewise  for  the  triple-excitation  correlation  energy. 

The  correlation  energy  from  the  single  excitations  within  subsystem  A  is  given  by 

£(fls)  =  2]TfV  (4.11) 

i&A 

where  the  summation  spans  the  occupied  orbitals  of  subsystem  A,  f'  is  the  internal- 
external  part  of  the  Fock  matrix  in  vector  form,  and  t*  are  the  single  excitation 
amplitudes  in  vector  form.47 

The  correlation  energy  from  the  double  excitations  within  subsystem  A  is  given 
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by 


A>)  =  £  (2  -  ^)tr[L«C«], 

i>j 

i,jeA 


(4.12) 


where  the  summation  spans  the  occupied  orbitals  of  subsystem  A,  and  L*-7  are  the 
internal  coulomb  and  exchange  matrices.  The  matrix  elements  of  are  given  by 
Cljs  =  T?s  +  frt{,  where  T^s  and  t'r  are  the  double  and  single  excitation  amplitudes, 
respectively.4' 

Finally,  the  correlation  energy  from  the  triple  excitations  within  subsystem  A  is 
given  by 

At)  =  £  (2  -  Sii  -  Sjt)  (  £  fr'Srr'  (j#<)  X% 

t>j>k  \rstr' 

i,j,k£A 

+  ^ Sss'  I kt )  Xrst  (4- 13^ 

rsts' 

+  £tfs„,  (ir\js)  x%  +  £  wgx% 

rstt'  rst 

where  the  first  summation  spans  the  occupied  orbitals  of  subsystem  A,  the  indices 
i,j,  k  represent  occupied  orbitals,  and  the  indices  r,  s,  t  represent  unoccupied  orbitals. 
Srri  is  an  element  of  the  overlap  matrix  of  the  projected  atomic  orbitals,  (ir\js)  are 
two-electron  integrals,  and  fr  are  the  single  amplitudes.  X ^  is  defined  as  Xlr^t  = 
4 T*Jst  —  2 —  2 Tl3Sr  —  2T*^  +  where  are  the  triples  amplitudes.  The 

tensor  element  contains  the  double-excitation  amplitudes,  50,51 


4.3.2  Calculation  Details 

All  geometry  optimizations  are  performed  using  Gaussian09  53  and  are  provided  in  the 
supplemental  information.  All  other  calculations  are  performed  in  Molpro  2012.1. 54 
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In  all  calculations  the  orbitals  are  localized  using  Pipek-Mezey  localization. 55  The 
atoms  associated  with  subsystem  A  for  each  reaction  are  given  in  the  supplemental 
information.  Any  localized  orbital  with  a  Lowdin  charge  of  0.4  on  an  atom  associated 
with  subsystem  A  is  included  in  the  set  of  orbitals  associated  with  subsystem  A.  All 
calculations  employ  a  level  shift  parameter  /i,  which  is  set  to  106  au.  All  KS-DFT 
calculations  employ  a  large  grid  for  the  exchange-correlation  functional  evaluation, 
achieved  by  specifying  the  Molpro  option  GRID =10” 10 .  For  computational  efficiency, 
all  LCCSD(T)  calculations  employ  density  fitting  (DF),  and  the  triples  are  approxi¬ 
mated  using  the  noniterative  (TO)  procedure. 50,51 

To  enable  the  rigorous  comparison  of  terms  from  the  LCCSD(T)  calculation  and 
the  embedding  calculation,  some  care  must  be  taken.  First,  all  orbital  pairs  are 
correlated  to  recover  the  energy  from  canonical  CCSD(T).  Second,  the  choice  of  or¬ 
bitals  must  be  consistent  between  the  LCCSD(T)  and  CCSD(T)-in-DFT  embedding 
calculations. 

In  the  WFT-in-DFT  embedding  method,  subsystem  B  comprises  KS  MOs,  and 
thus  evaluation  of  the  errors  resulting  from  using  the  DFT  energy  of  subsystem  B 
requires  the  use  of  KS  MOs  as  the  reference  MOs.  The  difference  between  canonical 
CCSD(T)  using  the  HF  reference  and  DF-LCCSD(TO)  using  the  KS  reference  is 
within  0.3  mAh  for  all  reactions  discussed  in  this  paper,  which  is  smaller  than  the 
errors  that  will  be  discussed  in  the  next  sections;  therefore,  in  the  next  section  we 
will  simply  refer  to  terms  calculated  from  DF-LCCSD(TO)  as  CCSD(T). 

Likewise,  consistent  evaluation  of  the  error  arising  from  the  embedding  potential 
requires  that  the  reference  MOs  for  the  embedded  CCSD(T)  calculation  on  subsystem 
A  be  obtained  from  the  corresponding  DFT  calculation.  This  selection  of  reference 
MOs  is  used  only  throughout  Sec.  Ill  C.  In  all  other  sections,  the  reference  MOs 
of  the  embedded  CCSD(T)  calculation  are  chosen  to  be  the  set  of  MOs  resulting 
from  an  embedded  HF  calculation.  We  note  that  the  difference  between  CCSD(T)- 
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in-DFT  embedding  where  the  MOs  for  subsystem  A  are  obtained  from  an  embedded 
DFT  calculation  compared  to  an  embedded  HF  calculation  is  within  0.2  rnP],  for  the 
reactions  considered  in  Sec.  Ill  C. 

Below  we  analyze  the  contributions  to  the  embedding  error  for  a  set  of  six  energies 
associated  with  different  reactions.  All  of  the  chosen  reactions  are  large  enough 
to  involve  partitioning  across  a  covalent  bond,  but  also  small  enough  to  allow  for 
calculation  of  the  full  CCSD(T)  reference  energy.  The  reactions  considered  are  given 

in  Table  I. 

The  data  set  consists  of  the  following  reactions:  (1)  activation  energy  for  the  sym¬ 
metric  Sn2  reaction  of  CP  and  propyl  chloride;  (2)  acid  hydrolysis  of  dimethylether 
to  form  methanol;  (3)  deprotonation  of  the  phenol  hydroxyl  group;  (4)  ring-closing 
isomerization  of  3-methylene-  1-heptene  to  form  butylcylobutane;  (5)  the  Diels- Alder 
reaction  of  2-methoxy- 1,3-butadiene  with  methyl  vinyl  ketone;  and  (6)  the  activation 
energy  for  the  Diels-Alder  reaction.  The  geometries  are  provided  in  the  supplemental 
information. 


reaction 

E  /  mAh 

1 

Sn2  activation  barrier 

7.8 

2 

acid  hydrolysis 

177.8 

3 

phenol  deprotonation 

568.8 

4 

ring  closing 

10.6 

5 

Diels-Alder  reaction 

63.1 

6 

Diels-Alder  barrier 

34.0 

Table  4.1:  CCSD(T)  reaction  energies  and  barriers  in  the  test  set  obtained  using 
cc-pVTZ  with  aug-cc-pV(T+d)Z  on  Cl,  and  aug-cc-pVTZ  for  all  atoms  for  reactions 
2-4. 56,57  For  ease  of  error  analysis,  we  adopt  a  sign  convention  in  which  all  reactions  or 
activation  processes  are  positive  in  energy.  Geometries  were  obtained  using  B3LYP 
with  6-311G*++  (reaction  1),  def2-TZVP  (reactions  2-4),  or  6-31G*  (reactions  5, 
6).58~61 


103 


4.3.3  Sources  of  Error  in  WFT-in-DFT  embedding 

4. 3. 3.1  Error  from  the  Embedding  Potential 
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Figure  4.1:  The  error  arising  from  the  embedding  potential  (bine  squares),  the  DFT 
energy  of  subsystem  B  (violet  triangles),  and  the  nonadditive  exchange-correlation 
energy  (green  diamonds)  compared  to  the  total  CCSD(T)-in-B3LYP  embedding  error 
(black  circles).  CCSD(T)  calculations  performed  on  the  full  system  are  used  as  the 
reference.  The  largest  source  of  error  is  the  nonadditive  exchange-correlation  energy 
functional. 

Now  we  discuss  how  comparison  of  terms  in  the  energy  expressions  for  CCSD(T) 
and  CCSD(T)-in-DFT  embedding  can  be  used  to  determine  the  error  arising  from 
the  embedding  potential.  The  energy  of  subsystem  A  from  the  CCSD(T)  calculation 
is  the  sum  of  the  HF  energy  (using  the  KS  density)  and  the  correlation  energy  of 
subsystem  A, 

-^CCSD(T)  =  ^hf[7A]  +  £fS)  +  -^(D)  +  ^(T)'  (4-14) 

The  total  energy  of  subsystem  A  from  a  CCSD(T)-in-DFT  embedding  calculation  is 

^emb  =  (4/A|lFA  in  B[7A,  7B]|TA) 

-tr[7A(hA“B[7\7B]-h)]. 


(4.15) 
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For  an  embedding  potential  that  includes  all  of  the  CCSD(T)  many-body  effects,  the 
energy  of  F'ccsd(t)  and  E^mh  would  be  identical;  therefore,  the  error  arising  from  the 
embedding  potential  is  calculated  as 


perror  jpA  77 'A 

-^pot  “  ^emb  -^CCSDCT)- 


(4.16) 


The  error  in  the  reaction  energies  arising  from  the  embedding  potential  is  therefore 
the  change  in  E^°T  between  products  and  reactants,  AE°T^°T. 

The  blue  squares  in  Figure  4.1  show  the  value  of  AE^*™  for  the  data  set,  compared 
to  the  total  CCSD(T)-in-B3LYP  embedding  error  shown  in  the  black  circles.  For  no 
system  is  the  error  larger  than  1.5  niFT ,  with  the  average  error  being  0.8  miy, .  This 
demonstrates  a  key  insight  of  this  paper,  which  is  that  the  embedding  potential 
calculated  using  WFT-in-DFT  embedding  is  very  accurate. 


4. 3. 3. 2  Error  from  Use  of  DFT  for  Subsystem  B 

Next,  we  quantify  the  WFT-in-DFT  embedding  error  resulting  from  treating  subsys¬ 
tem  B  using  DFT.  This  error  is  obtained  by  computing 


t^iB, error  _ 

-^DFT  — 


E^FT  [7B. 


—  (-Ehf[7B]  +  ^(B)  +  ^(D)  +  ^XT))  > 


(4.17) 


which  allows  for  a  direct  comparison  of  the  DFT  and  CCSD(T)  energies  of  subsystem 
B. 

The  values  calculated  for  AFjjp™1  are  shown  in  Figure  4.1  as  violet  triangles.  The 
largest  error  in  this  data  set  is  2.5  mf?h  and  the  average  error  is  1.5  mE^.  These  errors 
are  larger  than  those  resulting  from  the  embedding  potential,  but  are  still  relatively 
small  compared  to  the  total  WFT-in-DFT  embedding  error.  Therefore,  for  this  data 
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set,  DFT  does  an  adequate  job  describing  the  energy  change  localized  within  the 
environment  and  is  not  the  dominate  source  of  error. 

4. 3. 3. 3  Error  from  the  Nonadditive  Exchange-Correlation  Energy 

Finally,  we  analyze  the  error  that  arises  from  evaluation  of  the  nonadditive  exchange- 
correlation  energy  with  an  approximate  functional.  The  error  is  obtained  by  comput¬ 
ing 


=  BE,ft[7A.7B] 

(4.18) 

-  (£5?[7a,7b1  +  £S?D)  +  b^!<t))  . 

which  allows  for  the  direct  comparison  of  the  approximate  density  functional  to  the 
energy  obtained  at  the  CCSD(T)  level. 

The  values  for  AE°cd,error  are  given  in  Figure  4.1  as  green  diamonds.  This  term 
dominates  the  WFT-in-DFT  embedding  error,  with  the  largest  value  of  AE“(fd’error 
being  14.2  m£h,  and  the  average  value  being  7.2  miq, .  It  is  thus  this  term  that  is 
responsible  for  introducing  the  largest  error  in  the  WFT-in-DFT  embedding  method¬ 
ology. 

The  sum  of  AE®1™1’,  AEppp°r,  and  AE°(!ld’error  captures  all  of  the  discrepancy 
between  the  CCSD(T)-in-DFT  and  full  CCSD(T)  calculations.  Due  to  the  use  of 
density  fitting  and  the  noniterative  triples  approximation  used  in  the  CCSD(T)  cal¬ 
culation,  the  sum  of  these  errors  is  off  by  an  average  of  0.4  mEh  compared  to  the  total 
CCSD(T)-in-B3LYP  embedding  error;  this  makes  no  difference  in  the  interpretation 
of  the  data. 

To  confirm  that  our  results  are  not  sensitive  to  the  approximate  exchange-correlation 
functional,  we  repeated  the  analysis  using  both  PBE62  and  M0663  (not  shown).  These 
conclusions  are  robust  with  respect  to  the  approximate  exchange-correlation  func- 
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tional.  The  nonadditive  exchange-correlation  energy  remains  the  largest  source  of 
error,  followed  by  the  DFT  energy  of  subsystem  B.  Again,  we  find  that  DFT,  for  all 
of  the  functionals  tested,  provides  very  accurate  embedding  potentials. 

4.3.4  Improvement  of  the  Nonadditive  Exchange-Correlation 
Energy 


Figure  4.2:  Bar  graph  of  the  error  in  the  energy  obtained  from  CCSD(T)-in-B3LYP 
embedding  (black),  MP2-corrected  CCSD(T)-in-B3LYP  embedding  (red),  and  SOS- 
MP2-corrected  CCSD(T)-in-B3LYP  embedding.  CCSD(T)  calculations  performed  on 
the  full  system  are  used  as  the  reference. 

Having  determined  the  nonadditive  exchange-correlation  energy  to  be  the  dom¬ 
inate  source  of  error,  new  algorithms  can  be  proposed  to  calculate  this  term  more 
accurately.  One  approach  would  be  to  evaluate  the  nonadditive  exchange  exactly  and 
to  use  a  computationally  cheap  WFT  method,  such  as  MP2,64  to  evaluate  the  non¬ 
additive  correlation.  The  resulting  correction  to  the  WFT-in-DFT  embedding  energy 


is  then 
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EAUP2  =  En* i!-AFi  7B]  +  Y1  (2Trl  ~  TH)  K{ 

i&A  rs 
j£B 

(4.19) 

-^ft[7A,7B] 

-tr  [(7^F-7A)(hAinB[7A,7B]-h)], 

where  7^F  is  the  HF  embedded  density  of  subsystem  A,  T*i  is  the  MP2  amplitude,  and 
KlJs  are  the  exchange  two  electron  integrals.65  For  the  MP2  calculation  the  orbitals 
{c/7  }  U  {4>i}b  are  used,  which  allows  for  the  direct  calculation  of  the  MP2  correlation 
between  the  HF  orbitals  for  A  and  the  KS  orbitals  of  B. 

Figure  4.2  compares  the  CCSD(T)-in-B3LYP  embedding  error  (black)  to  the  MP2- 
corrected  CCSD(T)-in-B3LYP  embedding  error  (red).  The  average  error  of  WFT-in- 
DFT  embedding  is  4.6  mAh,  which  drops  to  1.2  mAh  when  the  MP2  correction  is 
applied.  Alternatively,  instead  of  calculating  the  full  MP2  energy  in  Eq.  4.19,  one 
could  only  calculate  the  scaled  opposite  spin  (SOS)-MP2  correlation  energy.66  Scaling 
the  opposite  spin  MP2  correlation  by  the  usual  empirical  factor  of  1.3  leads  to  the 
SOS-MP2-corrected  CCSD(T)-in-B3LYP  embedding  error  shown  in  blue  in  Figure 
4.2.  Applying  the  SOS-MP2  correction  results  in  an  average  error  of  1.1  mAh,  which 
is  essentially  the  same  error  as  that  of  the  full  MP2  correction,  and  only  requires 
computations  that  scale  N 4  compared  to  N5  for  the  full  MP2  energy. 

The  average  error  of  standard  MP2  calculations  on  these  systems  is  6.3  mAh  rel¬ 
ative  to  CCSD(T);  it  is  thus  clear  that  effectiveness  of  the  MP2  correction  does  not 
rely  on  the  MP2  energy  being  particularly  accurate  for  the  full  calculation.  Instead, 
we  observe  that  MP2  theory  accurately  represents  the  correlation  energy  between 
subsystems  A  and  B,  while  not  necessarily  representing  other  correlation  terms  accu¬ 
rately.  This  is  consistent  with  other  local  coupled-cluster  methods  that  treat  distant 
pairs  at  the  MP2  level. 4 ' 
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4.4  Results  II:  Continuity,  Convergence,  and  Con¬ 
jugation  in  WFT-in-DFT  embedding 

4.4.1  Potential  Energy  Surfaces 


1.0  1.5  2.0,  o  2,5  3.0 

Rc=c  /  A 


Figure  4.3:  (A)  Potential  energy  curves  for  the  dissociation  of  the  C-C  bond 

in  singlet  1-penten-l-one  obtained  using  CCSD(T)  (green),  KS-DFT  with  B3LYP 
(blue),  and  CCSD(T)-in-B3LYP  embedding  (black).  The  structure  was  reoptimized 
at  the  HF/cc-pVDZ  level  of  theory  for  each  value  of  the  C-C  bond  distance.56  The 
0=C=CH-  moiety  was  treated  at  the  CCSD(T)  level  for  the  CCSD(T)-in-B3LYP 
embedding  calculations.  (B-D)  The  error  in  CCSD(T)-in-B3LYP  embedding  (black) 
and  MP2-corrected  CCSD(T)-in-B3LYP  embedding  (red)  as  a  function  of  distance 
between  the  carbon-carbon  double  bond.  The  results  are  shown  for  three  partitionings 
of  the  molecule,  with  subsystem  A  corresponding  to  (B)  =C=CH-,  (C)  0=C=CH-, 
or  (D)  0=C=CH-CH2-. 

Next,  we  examine  the  potential  energy  surface  for  heterolytic  bond  cleavage.  Lo¬ 
cal  correlation  methods  show  discontinuities  in  the  potential  energy  surface  for  the 
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heterolytic  bond  cleavage  of  CO  dissociation  in  ketene.6'  Here,  we  study  a  related 
system,  CO  dissociation  in  1-penten-l-one. 

Panel  A  of  Figure  4.3  shows  potential  energy  curves  calculated  using  CCSD(T), 
B3LYP,  and  CCSD(T)-in-B3LYP  embedding.  The  cc-pVDZ  basis  was  used  for  all  cal¬ 
culations.  Here,  B3LYP  performs  relatively  well  near  equilibrium,  but  overestimates 
the  energy  by  up  to  16  miy,  near  dissociation.  The  CCSD(T)-in-B3LYP  calcula¬ 
tions  are  very  accurate  near  equilibrium  and  slightly  underestimate  the  energy  near 
the  dissociation  limit.  MP2-corrected  CCSD(T)-in-B3LYP  were  also  performed  for 
this  system;  the  results  are  not  shown  in  panel  A  of  Figure  4.3,  because  they  are 
graphically  indistinguishable  from  the  uncorrected  CCSD(T)-in-B3LYP  results. 

Panels  B  D  show  the  error  in  CCSD(T)-in-B3LYP  embedding  and  MP2-corrected 
CCSD(T)-in-B3LYP  embedding  for  three  different  subsystem  partitionings  of  the 
molecule.  The  error  and  the  change  of  the  slope  at  the  derivative  discontinuity  around 
1.5  A  decreases  by  treating  more  of  the  system  at  the  CCSD(T)  level.  Energy  dis¬ 
continuities  of  50  fiEh  are  seen  at  short  distances,  as  shown  in  Figure  4.7  of  the 
supplemental  information.  Like  other  local  correlation  methods,  abrupt  changes  in 
the  localized  orbitals  for  different  nuclear  configurations  lead  to  discontinuities  in  the 
WFT-in-DFT  embedding  energy  and  its  derivatives.  Here,  these  defects  are  small 
and  can  be  systematically  controlled  by  increasing  the  size  of  subsystem  A. 

4.4.2  WFT-in-DFT  Embedding  of  Conjugated  Systems 

A  demanding  case  for  any  embedding  methodology  is  the  partitioning  of  a  7r-conjugated 
system.  The  applicability  of  WFT-in-DFT  embedding  to  treat  such  systems  is  tested 
and  compared  to  systems  without  conjugation. 

First,  we  consider  the  dissociation  of  a  fluoride  anion  from  both  an  alkane  chain  (1- 
fluorodecane)  and  an  alkene  chain  (l-fluoro-l,3,5,7,9-decapentaene).  The  geometries 
for  both  compounds  and  their  dissociated  products  were  obtained  using  B3LYP/def2- 
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Figure  4.4:  (A)  The  error  in  CCSD(T)-in-B3LYP  embedding  (black  open  circles) 
and  MP2-corrected  CCSD(T)-in-B3LYP  embedding  (red  filled  circles)  as  a  function 
of  the  number  of  carbons  included  in  subsystem  A  for  the  dissociation  of  the  alkane. 
The  B3LYP  energy  is  given  by  the  black  dotted  line.  (B)  Contributions  to  the  WFT- 
in-DFT  error:  embedding  potential  (blue  open  circles),  DFT  for  subsystem  B  (violet 
filled  circles),  and  DFT  for  nonadditive  exchange-correlation  energy  (green  squares). 
(C)  DFT  Mulliken  population  of  the  density  associated  with  subsystem  B  on  the 
atoms  in  subsystem  A,  shown  for  1-fluorodecane  (black  open  circles)  and  the  dissoci¬ 
ated  alkane  chain  (red  filled  circles). 


TZVP.  All  CCSD(T)  and  embedding  calculations  were  performed  using  the  cc-pVDZ 
basis,  with  aug-cc-pVDZ  for  fluorine.56 

Figure  4.4A  shows  the  CCSD(T)-in-B3LYP  with  and  without  the  MP2  correction 
for  fluoride  anion  dissociation  from  the  alkane  chain.  Results  are  provided  for  a  num¬ 
ber  of  different  choices  of  the  subsystem  partitioning,  and  the  error  of  both  methods 
can  be  seen  to  rapidly  vanish  as  more  atoms  are  included  in  the  WFT  subsystem. 

The  individual  sources  of  error  in  the  CCSD(T)-in-B3LYP  embedding  calcula- 
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tions,  computed  in  the  same  way  as  in  Sec.  4.3.3,  are  shown  in  Figure  4.4B.  Again,  it 
is  observed  that  the  error  arising  from  the  embedding  potential  is  small,  accounting 
for  only  a  small  portion  of  the  total  error.  Unlike  previous  results,  the  error  arising 
from  treating  subsystem  B  at  the  DFT  level  is  of  similar  magnitude  as  the  nonadditive 
exchange-correlation  energy  error.  As  these  errors  are  of  opposite  sign,  evaluating  the 
nonadditive  exchange-correlation  energy  using  DFT  leads  to  a  favorable  cancellation 
of  error.  The  MP2  correction  only  increases  the  accuracy  of  the  subsystem  interac¬ 
tion  energy,  and  cannot  be  expected  to  correct  large  errors  associated  with  the  DFT 
energy  of  subsystem  B. 

Figure  4.4C  shows  the  Mulliken  population  of  the  density  associated  with  subsys¬ 
tem  B  on  the  atoms  associated  with  subsystem  A.  In  the  dissociated  product,  the 
density  associated  with  subsystem  B  distributes  onto  the  atoms  of  subsystem  A  to 
stabilize  the  positive  charge.  We  find  that  when  the  difference  of  this  quantity  is 
large  between  two  configurations,  there  is  typically  a  favorable  cancellation  of  error 
between  the  error  arising  from  treating  subsystem  B  using  DFT  and  the  error  arising 
from  evaluating  the  nonadditive  exchange-correlation  energy  using  DFT.  In  general, 
we  note  that  if  the  nonadditive  exchange-correlation  is  not  the  dominant  source  of 
error,  the  MP2  correction  cannot  significantly  improve  the  accuracy  of  the  embedding 
calculation. 

After  dissociation  of  the  fluoride  anion  from  l-fluoro-l,3,5,7,9-decapentaene,  the 
subsequent  geometry  optimization  leads  to  an  isomerization  where  the  proton  on  the 
second  carbon  moves  to  the  first.  Therefore,  the  analysis  for  this  reaction  begins  at 
the  second  carbon.  Figure  4.5A  shows  the  error  in  CCSD(T)-in-B3LYP  embedding 
(black  open  circles)  and  MP2-corrected  CCSD(T)-in-B3LYP  embedding  (red  filled 
circles)  as  a  function  of  the  number  of  carbons  included  in  subsystem  A  for  fluoride 
anion  dissociation  from  l-fluoro-l,3,5,7,9-decapentaene.  Unlike  the  alkane  case,  the 
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Figure  4.5:  (A)  The  error  in  CCSD(T)-in-B3LYP  embedding  (black  open  circles) 
and  MP2-corrected  CCSD(T)-in-B3LYP  embedding  (red  filled  circles)  as  a  function 
of  the  number  of  carbons  included  in  subsystem  A  for  the  dissociation  of  the  alkene. 
The  B3LYP  energy  is  given  by  the  black  dotted  line.  (B)  Contributions  to  the  WFT- 
in-DFT  error:  embedding  potential  (blue  open  circles),  use  of  DFT  for  subsystem  B 
(violet  filled  circles),  nonadditive  exchange-correlation  energy  (green  squares).  (C) 
DFT  Mulliken  population  of  the  density  associated  with  subsystem  B  on  the  atoms 
in  subsystem  A,  shown  for  l-fluoro-l,3,5,7,9-decapentaene  (black  open  circles)  and 
the  dissociated  alkene  chain  (red  filled  circles). 


alkene  case  exhibits  large  errors  which  slowly  decrease  once  the  majority  of  the  system 
is  treated  at  the  CCSD(T)  level. 

Figure  4.5B  shows  the  decomposition  of  the  contributions  to  the  error  in  CCSD(T)- 
in-B3LYP  embedding.  In  this  calculation,  the  error  arising  from  treating  subsystem 
B  using  DFT  is  the  dominate  source  of  error.  This  explains  why  the  error  remains 
large  until  the  majority  of  the  system  is  treated  at  the  CCSD(T)  level,  and  why  the 
MP2-correction  is  insufficient  to  reduce  the  error. 
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Figure  4.5C  shows  the  Mulliken  population  of  the  density  associated  with  subsys¬ 
tem  B  on  the  atoms  associated  with  subsystem  A  for  the  alkene  case.  As  with  the 
alkane  case,  a  large  difference  in  this  quantity  is  seen  between  the  flnorinated  and 
defluorinated  compounds.  This  observation  provides  insight  into  why  the  error  from 
the  DFT  energy  of  B  contributes  strongly  to  the  error  of  the  embedding  calculations. 

The  magnitude  of  the  change  in  the  dipole  moment  between  the  flnorinated  and 
defluorinated  compounds  is  show  in  Table  4.2  for  KS-DFT  with  B3LYP  and  CCSD. 
In  the  alkane  dissociation,  the  change  in  the  dipole  moment  is  large,  demonstrating 
a  small  polarizability,  and  there  is  good  agreement  between  KS-DFT  and  CCSD.  In 
the  alkene  dissociation,  the  change  in  dipole  moment  is  considerably  smaller  than 
the  alkane  case,  demonstrating  that  the  density  polarizes  to  stabilize  charge.  For 
the  alkene,  there  is  large  disagreement  between  KS-DFT  and  CCSD,  demonstrating 
the  known  failure  of  DFT  to  accurately  treat  polarizability  though  a  7r-conjugated 
system.68  Therefore,  when  there  are  large  errors  associated  with  KS-DFT,  these  large 
errors  will  affect  the  DFT  energy  of  subsystem  B,  causing  large  WFT-in-DFT  em¬ 
bedding  errors.  We  emphasize  that  for  cases  in  which  DFT  does  correctly  describe 
the  polarization  of  the  environment,  this  large  source  of  error  does  not  arise.  The 
failure  of  WFT-in-DFT  embedding  in  Figure  4.5  is  not  a  failure  of  embedding  itself, 
but  rather  a  failure  of  DFT  to  accurately  treat  the  polarizability  of  ^-conjugated 
systems. 


method 

dissociation 

exchange 

alkane 

B3LYP 

7.338 

0.781 

CCSD 

7.539 

0.802 

alkene 

B3LYP 

1.702 

0.551 

CCSD 

3.034 

0.630 

Table  4.2:  The  magnitude  of  the  change  in  the  dipole  moment  between  products 
and  reactants  for  the  dissociation  of  F~  from  the  alkane  and  alkene  chains,  as  well  as 
the  corresponding  magnitudes  for  the  H-F  exchange  reaction.  Values  are  reported  in 
atomic  units. 


114 


123456789 

Number  of  Carbons  in  Subsystem  A 


Figure  4.6:  The  error  in  CCSD(T)-in-B3LYP  embedding  (black  open  circles)  and 
MP2-corrected  CCSD(T)-in-B3LYP  embedding  (red  filled  circles)  as  a  function  of  the 
number  of  carbons  included  in  subsystem  A  for  the  exchange  of  fluoride  to  a  hydride 
in  (a)  1-fluorodecane  and  (b)  l-fluoro-l,3,5,7,9-decapentaene. 


Finally,  we  consider  the  reaction  of  exchanging  the  fluoride  anion  from  1-fluorodecane 
and  l-fluoro-l,3,5,7,9-decapentaene  with  a  hydride  (Figure  4.6).  The  change  in  dipole 
moment  for  these  reactions  is  provided  in  Table  4.2.  These  reactions  exhibit  a  mod¬ 
erate  change  in  dipole  moment,  and  there  is  good  agreement  between  CCSD  and 
KS-DFT. 

Figures  4.6A  and  4.6B  plot  the  error  in  the  CCSD(T)-in-B3LYP  embedding  and 
MP2-corrected  CCSD(T)-in-B3LYP  embedding  energies  for  the  hydride  exchange  re¬ 
actions  from  alkane  and  alkene  chains,  respectively,  as  a  function  of  the  number  of 
carbons  included  in  subsystem  A.  For  every  partition,  the  errors  are  small.  For  the 
smallest  division,  the  MP2  correction  provides  a  significant  improvement  in  the  accu- 
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racy  of  the  CCSD(T)-in-B3LYP  embedding  energy;  for  larger  divisions,  the  effect  of 
the  MP2  correction  is  much  smaller.  Unlike  in  the  case  of  fluoride  anion  dissociation, 
DFT  applied  to  the  hydride  exchange  reaction  accurately  represents  the  change  in 
dipole.  As  there  are  no  large  errors  arising  from  the  DFT  energy  of  subsystem  B, 
WFT-in-DFT  embedding  performs  accurately  and  the  MP2  correction  further  im¬ 
proves  the  energetics. 

The  important  observation  from  these  calculations  is  that  when  there  is  a  large 
error  in  the  DFT  calculation  on  the  environment,  there  will  be  correspondingly  large 
errors  in  the  WFT-in-DFT  embedding  energy.  Importantly,  this  failure  is  associated 
with  errors  intrinsic  to  the  DFT  functionals,  and  does  not  arise  due  to  errors  in  the 
embedding  potential.  When  a  chemical  process  involves  a  large  change  in  the  Mulliken 
population  of  subsystem  B  located  on  the  subsystem  A  atoms,  it  is  likely  that  the 
embedding  error  will  be  dominated  by  errors  arising  from  the  DFT-level  treatment 
of  subsystem  B;  errors  of  this  sort  cannot  be  reduced  by  the  MP2  correction. 

4.5  Conclusions 

Projector-based  quantum  embedding  provides  a  scheme  for  multiscale  descriptions 
with  the  exactness  property  that  DFT-in-DFT  is  equivalent  to  DFT  on  the  whole 
system.36,3'  In  many  tests  and  applications,  we  find  the  accuracy  of  the  scheme  to 
be  excellent,  allowing  for  aggressive  partitioning  across  covalent  bonds  close  to  the 
reactive  center  of  the  system  of  interest.  However,  for  some  applications,  the  errors 
introduced  by  embedding  are  larger  than  would  typically  be  acceptable,  and  the 
principal  aims  of  this  paper  have  been  to  understand  and  take  steps  towards  resolving 
the  errors  in  such  cases. 

Careful  comparison  of  CCSD(T)-in-DFT  embedding  calculations  with  full  CCSD(T) 
calculations  has  led  to  key  insights  regarding  the  sources  of  error  in  the  embedding 
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calculations.  First,  the  embedding  potential  obtained  using  approximate  density  func¬ 
tionals  is  found  to  be  accurate  for  all  of  the  cases  we  have  investigated,  making  a 
contribution  to  the  overall  error  of  the  embedding  calculation  that  is  negligible  com¬ 
pared  to  other  sources  of  error.  It  was  not  immediately  obvious  that  this  would  be 
the  case,  because  functionals  (particularly  in  cases  where  they  are  parameterized)  are 
designed  with  accurate  energies  in  mind. 

And  second,  it  is  found  that  in  many  cases,  the  primary  source  of  error  in 
CCSD(T)-in-DFT  embedding  is  the  treatment  of  nonadditive  exchange-correlation 
effects  with  an  approximate  density  functional.  This  is  important  because  it  is  the 
one  term  in  the  error  for  which  simple  corrections  can  be  developed  that  conserve  the 
efficiency  of  the  original  method.  Here,  we  found  that  use  of  MP2  or  SOS-MP2  cor¬ 
rections  for  this  term  typically  improved  the  accuracy  of  the  energetics  for  chemical 
reactions,  reducing  the  average  error  from  4.6  mAh  to  1.2  mAh  with  respect  to  full 
CCSD(T)  calculations. 

To  investigate  the  convergence  with  respect  to  the  size  of  subsystem  A,  we  studied 
dissociation  and  exchange  events  at  the  terminus  of  10-carbon  alkyl  and  conjugated 
chains.  For  the  removal  of  F~,  the  results  of  the  CCSD(T)-in-DFT  embedding  cal¬ 
culation  for  the  conjugated  system  are  noticeably  worse  than  for  the  alkane,  and  it  is 
found  that  the  MP2  correction  does  not  reduce  this  error  in  the  computed  reaction 
energy.  Our  analysis  shows,  however,  that  these  results  follow  from  the  fact  that 
DFT  provides  a  poor  description  of  the  polarization  of  the  charged  alkene  fragment 
and  that  the  uncorrected  CCSD(T)-in-DFT  results  benefit  from  a  cancellation  of  er¬ 
rors  in  the  DFT  treatment  of  subsystem  B  and  in  the  DFT  treatment  of  nonadditive 
exchange-correlation.  The  MP2  correction  improves  the  description  of  nonadditive 
energy  term,  but  it  does  not  compensate  for  the  inaccuracies  in  the  DFT  description 
of  subsystem  B. 

For  a  hydride  exchange  reaction  at  the  terminus  of  the  alkyl  and  conjugated  chains, 
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the  CCSD(T)-in-DFT  embedding  results  converge  smoothly  and  rapidly  to  reference 
CCSD(T)  calculations  performed  over  the  full  system,  regardless  of  inclusion  of  the 
MP2  correction  and  regardless  of  conjugation  in  the  chain.  These  results  demonstrate 
that  in  the  regime  where  DFT  is  adequate  for  the  treatment  of  the  environment, 
our  projector-based  embedding  scheme  can  effectively  partition  the  system,  even  in 
conjugated  molecules. 

The  current  work  demonstrates  that  projection-based  embedding  provides  both 
a  rigorous  and  practical  approach  to  embedding  correlated  wavefunctions  in  a  DFT 
description  of  the  environment.  Although  the  results  presented  here  utilize  coupled- 
cluster  methods  for  describing  the  correlated  wavefunction,  we  emphasize  that  projection- 
based  embedding  can  be  combined  just  as  easily  with  multi-reference  electronic  struc¬ 
ture  methods,  as  well  as  any  mean-field  description  of  the  environment.  The  embed¬ 
ding  method  is  straightforward  to  employ  -  requiring  only  the  specification  of  which 
atoms  are  to  be  treated  at  the  WFT  and  DFT  levels  of  theory  -  and  it  is  fully  imple¬ 
mented  and  available  in  the  Molpro  quantum  chemistry  package. 

4.6  Accurate  and  Systematically  Improvable  Den¬ 
sity  Functional  Theory  Embedding  for  Corre¬ 
lated  Wavefunctions:  Supplemental  Informa¬ 
tion 

4.6.1  Potential  Energy  Surfaces 
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Figure  4.7:  Graph  of  the  error  in  CCSD(T)-in-B3LYP  (black)  and  MP2  corrected 
CCSD(T)-in-B3LYP  (red)  as  a  function  of  distance  between  the  carbon-carbon  dou¬ 
ble  bond.  The  0=C=CH-CH2-  moiety  was  treated  at  the  CCSD(T)  level  for  the 
CCSD(T)-in-B3LYP  embedding  calculations.  Abrupt  changes  in  the  localized  or¬ 
bitals  for  different  nuclear  configurations  lead  to  discontinuities  in  the  WFT-in-DFT 
energy  and  its  derivatives 


4.7  Data  Set  Computational  Details 


Molecular  geometries  are  reported  in  angstrom  using  Cartesian  coordinates. 

4.7.1  Symmetric  Sn2  Reaction  Barrier 


Propylchloride 

KS-DFT  (B3LYP)  /  6-311g*++  Optimized  Geometry 
Subsystem  A  atoms:  C8,  H9,  H10,  Clll 
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Cl 

-1.823996 

-0.717480 

0.144292 

H2 

-2.784508 

-0.872856 

-0.353495 

H3 

-1.203709 

-1.595688 

-0.046137 

H4 

-2.013376 

-0.669699 

1.221397 

C5 

-1.151988 

0.559917 

-0.359191 

H6 

-1.810375 

1.420026 

-0.177162 

H7 

-1.004107 

0.511737 

-1.442695 

C8 

0.173632 

0.884185 

0.310894 

H9 

0.573039 

1.838983 

-0.026992 

H10 

0.089151 

0.898664 

1.397499 

Clll 

1.468705 

-0.346524 

-0.067552 

Transition  State 


KS-DFT  (B3LYP)  /  6-311g*++  Optimized  Geometry 
Subsystem  A  atoms:  Cl,  H2,  H3,  Clll,  012 


Cl 

0.000719 

-0.388805 

0.105519 

H2 

0.000666 

-1.250604 

0.746120 

H3 

0.000734 

-0.559306 

-0.953281 

C4 

0.000193 

0.996697 

0.687217 

H5 

-0.878266 

1.089608 

1.331323 

H6 

0.879535 

1.090788 

1.329810 

C7 

-0.001812 

2.094992 

-0.373882 

H8 

-0.005528 

3.088594 

0.087416 

H9 

0.887844 

2.018883 

-1.005190 

H10 

-0.889460 

2.013099 

-1.006979 

Clll 

2.416740 

-0.696767 

-0.089570 

C112 

-2.416159 

-0.697843 

-0.089391 
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4.7.2  Acid  Hydrolysis  Reaction  Energy 


Dimethyl  Ether 

KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 
Subsystem  A  atoms:  Cl,  02,  H3,  H4,  H5 


Cl 

-1.173920 

0.194451 

0.000000 

02 

-0.000001 

-0.586419 

0.000000 

H3 

-1.231875 

0.834561 

-0.890754 

H4 

-1.231865 

0.834609 

0.890722 

H5 

-2.021089 

-0.490198 

0.000022 

C6 

1.173921 

0.194451 

0.000000 

H7 

2.021088 

-0.490200 

0.000005 

H8 

1.231861 

0.834579 

0.890742 

H9 

1.231879 

0.834589 

-0.890736 

Methanol 

KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 
Subsystem  A  atoms:  02,  H6 


Cl 

0.665448 

-0.020220 

0.000000 

02 

-0.748210 

0.122194 

0.000000 

H3 

1.028648 

-0.544325 

-0.891167 

H4 

1.028648 

-0.544324 

0.891167 

H5 

1.083163 

0.985675 

0.000000 

H6 

-1.147468 

-0.753259 

0.000000 

ch3+ 

KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 

Subsystem  A  atoms:  Cl,  H2, 

H3,  H4 

Cl 

0.000000 

-0.000001 

0.000000 

H2 

0.945569 

-0.545889 

-0.000001 

H3 

-0.945544 

-0.545931 

-0.000001 

H4 

-0.000024 

1.091829 

-0.000001 
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4.7.3  Phenol  Deprotonation 

Phenol 


KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 
Subsystem  A  atoms:  01,  H6 


Cl 

0.218552 

-1.219283 

0.000014 

C2 

-1.168651 

-1.184791 

-0.000069 

C3 

-1.850797 

0.028852 

-0.000109 

C4 

-1.127187 

1.214896 

-0.000067 

C5 

0.263118 

1.193734 

0.000013 

C6 

0.936484 

-0.025561 

0.000056 

HI 

0.758987 

-2.15667 

0.00005 

H2 

-1.721496 

-2.115963 

-0.000101 

H3 

-2.932538 

0.048567 

-0.000172 

H4 

-1.643317 

2.16685 

-0.000098 

H5 

0.823932 

2.12249 

0.00004 

01 

2.300596 

-0.110893 

0.000128 

H6 

2.680547 

0.774789 

0.000232 

Phenolate 


KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 
Subsystem  A  atoms:  01 


Cl 

1.077309 

0.000000 

-0.000014 

C2 

0.286113 

1.209101 

-0.000003 

C3 

-1.097619 

1.196760 

0.000012 

C4 

-1.823396 

0.000000 

0.000013 

C5 

-1.097619 

-1.196759 

0.000010 

C6 

0.286113 

-1.209102 

-0.000005 

HI 

0.829304 

2.149743 

-0.000005 

H2 

-1.634545 

2.143825 

0.000017 

H3 

-2.907447 

0.000000 

0.000035 

H4 

-1.634545 

-2.143825 

0.000013 

H5 

0.829303 

-2.149744 

-0.000009 

01 

2.341566 

0.000000 

-0.000016 
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4.7.4  Ring  Closing 


Ring  Open 

KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 
Subsystem  A  atoms:  012,015, C16,H17,H18>C19,H20,H21,H22 


Cl 

-3.545023 

0.629867 

-0.421139 

02 

-2.506770 

-0.227213 

0.300475 

H3 

-3.361405 

1.694576 

-0.256667 

H4 

-4.555825 

0.409794 

-0.071888 

H5 

-3.521013 

0.454536 

-1.499548 

06 

-1.075186 

0.045717 

-0.160699 

H7 

-2.739265 

-1.286652 

0.148181 

H8 

-2.577664 

-0.054295 

1.379802 

09 

-0.033543 

-0.824403 

0.557100 

H10 

-0.838343 

1.103391 

-0.006521 

Hll 

-0.998423 

-0.127598 

-1.239341 

C12 

1.388709 

-0.552578 

0.125819 

H13 

-0.266039 

-1.879491 

0.392857 

H14 

-0.114525 

-0.651038 

1.637376 

C15 

1.939021 

0.764969 

0.494120 

C16 

2.110936 

-1.464774 

-0.530108 

H17 

3.144393 

-1.286877 

-0.798480 

H18 

1.691125 

-2.425795 

-0.801128 

C19 

2.771005 

1.502943 

-0.237972 

H20 

1.601127 

1.161132 

1.449678 

H21 

3.142039 

2.454299 

0.121039 

H22 

3.098937 

1.186855 

-1.220934 
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Ring  Closed 

KS-DFT  (B3LYP)  /  Def2-TZVP  Optimized  Geometry 
Subsystem  A  atoms:  C12,C15,C16,H17, H18,C19, H20,H21,H22 


Cl 

3.794291 

0.519261 

0.000003 

C2 

2.639877 

-0.480728 

-0.000003 

H3 

3.758022 

1.164438 

-0.881470 

H4 

4.761328 

0.012414 

-0.000002 

H5 

3.758024 

1.164425 

0.881486 

C6 

1.263955 

0.185711 

0.000001 

H7 

2.724185 

-1.134284 

0.874923 

H8 

2.724184 

-1.134273 

-0.874938 

C9 

0.107415 

-0.813847 

-0.000002 

H10 

1.176665 

0.837654 

-0.875560 

Hll 

1.176667 

0.837647 

0.875568 

C12 

-1.248215 

-0.197709 

-0.000001 

H13 

0.192433 

-1.477758 

0.870558 

H14 

0.192433 

-1.477754 

-0.870565 

C15 

-1.733584 

1.050750 

-0.000003 

C16 

-2.624861 

-0.841354 

0.000003 

H17 

-2.872972 

-1.429596 

-0.887456 

H18 

-2.872969 

-1.429592 

0.887466 

C19 

-3.185495 

0.621333 

0.000000 

H20 

-1.268168 

2.028627 

-0.000006 

H21 

-3.765065 

0.888772 

-0.886974 

H22 

-3.765062 

0.888776 

0.886975 
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4.7.5  Diels  Alder 


Methyl  Vinyl  Ketone 

KS-DFT  (B3LYP)  /  6-31G*  Optimized  Geometry 
Subsystem  A  atoms:  C1,H1>C2,H2,H3,C3,01 


Cl 

0.876535 

-0.634977 

-0.000004 

HI 

0.993976 

-1.716638 

0.000005 

C2 

1.942550 

0.172493 

0.000004 

H2 

2.953994 

-0.224441 

0.000020 

H3 

1.852200 

1.255508 

-0.000004 

C3 

-0.544547 

-0.189222 

-0.000028 

01 

-1.435445 

-1.025231 

0.000011 

C4 

-0.858775 

1.298350 

0.000004 

H4 

-0.434414 

1.789750 

0.883609 

H5 

-0.434444 

1.789780 

-0.883599 

H6 

-1.942326 

1.428021 

0.000024 

2-  met  hoxy- 1 , 3-butadiene 

KS-DFT  (B3LYP)  /  6-31G*  Optimized  Geometry 

Subsystem  A  atoms: 

C1,C4,C2,H3,C3,H4,H5,H1jH2 

Cl 

0.617061 

1.646444 

-0.000018 

HI 

1.694859 

1.746315 

-0.000110 

H2 

0.038603 

2.562692 

0.000138 

C2 

-1.487945 

0.367430 

0.000047 

C3 

-2.185049 

-0.774091 

0.000042 

H3 

-2.005313 

1.324367 

0.000116 

H4 

-3.270785 

-0.767429 

0.000099 

H5 

-1.692300 

-1.740640 

-0.000031 

01 

0.579814 

-0.762164 

-0.000208 

C4 

-0.023554 

0.461802 

-0.000068 

C5 

1.997125 

-0.795343 

0.000129 

H6 

2.407626 

-0.307484 

-0.893910 

H7 

2.407176 

-0.307121 

0.894173 

H8 

2.275794 

-1.850840 

0.000397 
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Transition  State 

KS-DFT  (B3LYP)  /  6-3 1G*  Optimized  Geometry 
Subsystem  A  atoms:  C1,C2,C3,C4,C5,C6,H1)H2,H3,H4,H5,H6,H7,H8>C7,01 


Cl 

1.068810 

-0.812570 

0.719810 

C2 

-0.339230 

-1.188850 

-0.695950 

C3 

-1.266210 

-0.148460 

-0.876330 

C4 

-0.304900 

1.715430 

0.605970 

C5 

1.004030 

1.590610 

0.231600 

C6 

1.712240 

0.368160 

0.321740 

HI 

-1.130060 

0.556110 

-1.689640 

H2 

0.355970 

-1.381110 

-1.508990 

H3 

-0.678730 

-2.097060 

-0.202640 

H4 

1.635290 

-1.737950 

0.762130 

H5 

-0.866540 

2.614980 

0.375910 

H6 

-0.764700 

1.054020 

1.327930 

H7 

1.487360 

2.384380 

-0.333520 

H8 

0.347890 

-0.721890 

1.524310 

C7 

-2.592740 

-0.095470 

-0.255940 

01 

-3.435810 

0.724000 

-0.615230 

C8 

-2.932550 

-1.097300 

0.849810 

H9 

-3.061010 

-2.104130 

0.432560 

H10 

-2.148450 

-1.163030 

1.614860 

Hll 

-3.870750 

-0.793220 

1.318090 

02 

2.936890 

0.401990 

-0.277740 

C9 

3.765940 

-0.749750 

-0.211570 

H12 

4.727800 

-0.448890 

-0.630610 

H13 

3.904990 

-1.083200 

0.824760 

H14 

3.359950 

-1.577760 

-0.806160 
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Product 

KS-DFT  (B3LYP)  /  6-3 1G*  Optimized  Geometry 
Subsystem  A  atoms:  C3,C8,C2,C1,C4,C5,H7,H6,H8>H5,H2,H3,H1>C6,01 


Cl 

-1.097853 

-0.009400 

0.307622 

HI 

-0.934291 

-0.461527 

1.296642 

C2 

-0.395447 

-0.905157 

-0.741185 

H2 

-0.861748 

-1.897406 

-0.785253 

H3 

-0.527828 

-0.444533 

-1.727887 

C3 

1.101813 

-1.049944 

-0.432388 

H4 

1.246067 

-1.766230 

0.390560 

H5 

1.616769 

-1.478343 

-1.302945 

C4 

-0.485038 

1.400083 

0.279829 

H6 

-0.903187 

1.945611 

-0.577138 

H7 

-0.808110 

1.957516 

1.170221 

C5 

1.019771 

1.371473 

0.219750 

H8 

1.558008 

2.294164 

0.423611 

C6 

-2.601257 

0.012621 

0.023153 

01 

-3.092155 

0.811145 

-0.753610 

C7 

-3.449014 

-1.027538 

0.734995 

H9 

-3.495173 

-0.790941 

1.806419 

H10 

-4.460915 

-1.032016 

0.324814 

Hll 

-3.005536 

-2.027088 

0.650795 

02 

3.102253 

0.393352 

-0.206064 

C8 

1.730467 

0.282418 

-0.094316 

C9 

3.875760 

-0.623808 

0.415869 

H12 

3.702491 

-1.613323 

-0.027993 

H13 

4.921497 

-0.349057 

0.259886 

H14 

3.675958 

-0.677301 

1.495684 
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Outlook 

This  work  shows  the  significant  progress  contributed  by  this  thesis  to  the  development 
of  embedded  density  functional  theory  for  the  treatment  of  correlated  wavefunctions. 
We  started  with  a  simple  single  electron  embedded  in  an  atomic  system,  progressed 
to  covalent  bonded  molecular  systems,  and  advanced  to  transition  metals  and  large 
complex  chemical  systems.  The  methodology  has  not  only  advanced  the  size  and 
complexity  of  systems  which  can  be  treated,  but  also  the  accuracy  and  robustness. 
However,  significant  challenges  remain  to  be  addressed  and  are  discussed  briefly  here. 

Chapter  4  discussed  the  use  of  computationally  inexpensive  wavefunction  methods 
to  evaluate  the  nonadditive  exchange-correlation  energy.  In  those  cases  we  applied  the 
use  of  MP2  theory,  and  it  proved  very  accurate.  However,  there  are  many  systems  in 
which  MP2  theory  will  give  qualitatively  inaccurate  results,  such  as  transition  metal 
complexes  and  complexes  with  near  degeneracies.  Therefore,  alternative  inexpensive 
wavefunction  methods,  such  as  the  random  phase  approximation,  should  be  explored 
to  determine  the  accuracy  of  those  methods  in  the  context  of  WFT-in-DFT. 

The  methodology  used  in  chapter  4  employs  the  KS-DFT  density  for  subsystem 
B,  which  was  obtained  from  performing  a  KS-DFT  calculation  on  the  full  system. 
However,  in  some  systems  the  density  in  subsystem  A  may  change  significantly  from 
KS-DFT  to  wavefunction  theory.  Therefore,  the  next  challenge  to  the  field  is  how  to 
relax  the  density  of  subsystem  B,  after  one  obtains  the  density  of  subsystem  A.  This 
self-consistent  minimization  allows  for  the  two  subsystems  to  be  correctly  polarized. 

This  self-consistent  procedure  will  allow  for  the  total  systems  orbitals  to  be  at  the 
absolute  minimum,  which  has  several  advantages.  First,  it  allows  for  the  treatment  of 
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excited  states.  If  an  excited-state  wavefunction  calculation  is  performed  on  subsystem 
A,  subsystem  B  can  polarize  to  the  excited  state  and  should  be  significantly  more 
accurate  than  frozen  density  embedding. 

A  self-consistent  procedure  would  also  allow  for  the  calculation  of  atomic  forces 
through  the  use  of  the  Hellmann-Feynman  theorem.  This  would  allow  for  geometry 
optimizations  at  the  WFT-in-DFT  level  of  theory.  Previous  experience  has  shown 
that  geometries  obtained  at  the  KS-DFT  level  of  theory  can  be  inaccurate  compared 
to  experimental  geometries,  particularly  for  transition  metal  complexes.  However, 
this  requires  the  use  of  WFT  methods  with  computationally  tractable  gradients.  For 
many  WFT  methods,  geometry  optimization  is  too  costly  for  all  but  the  smallest  of 
systems,  and  this  remains  a  great  challenge  to  the  WFT  community. 

Although  several  challenges  remain  in  the  WFT-in-DFT  methodology,  it  is  im¬ 
portant  to  emphasize  how  breakthroughs  in  both  WFT  and  DFT  can  be  seamlessly 
and  effortlessly  integrated  into  the  held.  In  this  sense,  this  work  will  remain  compli¬ 
mentary  to  the  work  being  done  in  both  the  WFT  and  DFT  scientific  communities. 


